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1.  Introduction 

With  the  term  "gravity  field  modelling"  we  usually,  in  the  geodetic  community 
mean  methods  for  representing  the  external  potential  of  the  earth,  in  order  to 
be  able  to  estimate  quantities  related  to  the  gravity  field  from  a  given  set  of 
"observed"  quantities.  Such  methods  include  spherical  harmonic  expansions,  integral 
formulas  such  as  Stoke;'  and  Vening-Meinesz '  formulas  and  "spatial"  approximation 
methods  such  as  collocition  and  generalized  point  mass  modelling  (e.g.  Bjerhammar's 
methods).  Common  for  ill  of  these  methods  is  that  they  are  approximation  methods 
for  harmonic  functions,  all  rely  on  the  assumption  that  the  anomalous  potential 
T  fulfills  Laplace's  equation  V  2T  =  0  at  least  outside  the  surface  of  the  earth. 

No  assumption  is  made  regarding  the  density  distribution  actually  generating  the 
gravity  field. 

In  contrast  the  term  "gravity  field  modelling"  as  used  in  geophysics  stands 
for  the  process  of  determining  internal  density  distributions  of  the  earth,  consis¬ 
tent  with  the  observed  outer  field.  This  inversion  is  non-unique,  and  to  get 
reasonable  answers  the  geophysicist  must  introduce  constraints,  through  selection 
of  a  finite  dimensional  representation  of  the  structures  (density  values,  depth 
parameters,  etc.),  through  "fixing"  some  of  these  parameters  based  on  independent 
geophysical  information  (well  data,  seismic  interpretations)  and  through  subjective 
choices  of  the  most  "realistic"  models  in  terms  of  geology.  At  the  basis  of  the 
geophysical  gravity  field  modelling  is  the  "direct  problem"  of  potential  field 
theory:  to  calculate  the  gravity  potential  and  its  derivatives  in  space  due  to 
given  density  distributions. 

When  the  prime  interest  is  in  "external"  gravity  field  modelling,  any  geophysical 
density  model,  realistic  or  not,  may  in  principle  be  used  to  represent  a  part 
of  the  external  field  through  a  direct  computation  of  the  effects  of  the  assumed 
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density  distribution.  If  the  density  distribution  is  realistic  we  would  expect 
the  remaining  field  to  be  more  smooth,  in  some  cases  the  fit  would  be  so  close 
that  we  would  have  no  use  for  any  external  gravity  field  modelling  at  all.  Usually, 
however,  our  knowledge  of  density  anomalies  is  confined  to  more  shallow  structures 
of  crustal  and  upper  mantle  origin,  thus  mainly  contributing  to  the  shorter  wavelengths 
of  the  variation  of  the  gravity  field.  Longer  wavelength  parts  of  the  signal  are 
more  conveniently  treated  using  "external"  modelling,  such  as  high  degree  and 
order  spherical  harmonic  expansions  of  the  geopotential. 

The  "external"  and  "internal"  modelling  may  conveniently  be  done  simultaneously 
in  some  cases,  using  e.g.  combined  versions  of  collocation  and  geophysical  inversion 
procedures.  Thus  we  will  at  the  same  time  estimate  both  the  external  gravity 
field  and  density  values  inside  the  earth.  This  approach  has  the  advantage  that 
it  allows  fairly  easy  use  of  independent  geologic/geophysical  information  as  data 
for  the  construction  of  external  gravity  field  models.  Due  to  the  difficult  choice 
of  geologically  "reasonable"  density  parameters  it  will,  however,  hardly  ever 
be  a  "hands-off"  automatic  method  like  standard  collocation. 

Combined  collocation/inversion  will  probably  prove  itself  useful  for  inversion 
of  future  "multi sensor"  gravity  data,  as  e.g.  gravity  vector  measurements  by  • 
inertial  survey  systems  and  gravity  gradiometer  measurements.  In  conventional 
geophysical  inversion  we  have  an  inherent  arbitrary  choice  of  the  "regional"  field, 
representing  the  effects  of  all  other  sources  than  the  density  structures  of  interest.. 
This  regional/residual  -  separation  is  done  using  more  or  less  crude  forms  of 
filtering  and  trend  fitting.  When  we  have  several  different  types  of  gravity 
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field  data  containing  information  about  the  same  mass  body,  it  is  essential  that 
the  filtering  is  consistent  between  the  various  gravity  field  data  types,  such 
that  the  filter  output  still  represent  quantities  related  to  the  same  harmonic 
function.  This  will  be  assured  using  combined  collocation/inversion  and  similar 


methods. 


In  this  report  tie  utilization  of  density  anomalies  for  genereral ized  gravity 


field  modelling  will  be  treated  in  the  first  chapters  in  a  rather  broad  way. 

The  bulk  of  the  report  will,  however,  be  concerned  about  the  most  important  and 
also  best  known  density  anomalies  on  the  earth,  namely  density  anomalies  associated 
with  the  topography. 

The  density  anomalies  relating  to  the  topography  include  the  direct  gravita¬ 
tional  effects  of  the  visual  topography  on  the  continents,  the  ocean  bathymetry, 
ice  cap  effects  and  the  isostatic  compensation.  These  effects  together  represent 
a  major  part  of  the  variation  of  the  earth's  gravity  field,  especially  at  shorter 
wavelengths  (less  than,  say,  a  few  hundred  kilometers),  where  the  direct  computed 
topographic  effects  only  to  a  low  degree  are  counteracted  by  the  isostatic  compensa 
tion  effects.  In  mountainous  areas  the  topographic  effects  completely  dominate 
the  local  variation  of  the  gravity  field,  and  some  kind  of  terrain  reduction  is 
indispensable  when  attempting  gravity  field  modelling  in  such  areas.  The  most 
well-known  terrain  reduction  is  the  Bouguer  reduction,  which  is  well-suited  for 
geophysical  work  and  prediction  of  mean  free-air  anomalies  (for  use  e.g.  in  tradi¬ 
tional  geodetic  gravity  field  modelling),  but  is  not  applicable  for  reduction 
of  geoid  undulations.  Isostatic  reductions  provide  the  smoothest  possible  residual 
fields  on  a  global  basis,  and  are  easily  applicable  to  all  the  various  types  of 
gravity  field  data  available.  The  computation  of  topographic/isostatic  effects 
is  facilitated  by  high-degree  spherical  harmonic  expansions  of  the  isostatic  reduc¬ 
tion  potential  (Rapp,  1982),  but  for  local  applications  the  computations  are  still 
relatively  laborious.  Since  the  usual  Airy-type  isostasy  does  not  operate  on 
a  local  scale  (the  short  wavelength  topography  is  supported  by  the  finite  strength 
of  the  earth's  crust)  we  might  not  like  to  introduce  the  somewhat  arbitrary  isosta 
tic  reduction  mass  anomalies  at  the  crust/mantle  boundary.  Instead  we  might  just 
simply  try  to  take  only  the  short  wavelength  variations  of  the  topography  into 


account, introducing  an  arbitrary  smooth  mean  elevation  surface  ("model"  earth) 
as  a  "reference"  topography,  the  gravitational  influence  of  which  is  not  directly 
taken  into  account  in  the  gravity  field  modelling.  For  gravity  anomalies  such 
a  residual  terrain  correction  corresponds  closely  to  the  formation  of  regional 
mean  free-air  anomalies,  and  by  choice  of  a  proper  reference  elevation  surface, 
such  as  defined  through  a  spherical  harmonic  expansion  of  the  earth's  topography 
to  degree  and  order  180,  we  end  up  with  a  reduction  which  would  be  expected  to 
give  somewhat  similar  results  as  conventional  isostatic  reductions. 

For  local  modelling  of  the  gravity  field  -  on  which  the  main  emphasis  is 
put  in  this  report  -  the  availability  of  high  degree  and  order  spherical  harmonic 
expansions  of  the  geopotential  (Lerch  et  al.,  1981;  Rapp,  1982)  has  proven  itself 
to  be  a  major  break-through  of  big  practical  value.  For  a  region  like  Scandinavia 
with  reliable  l°x  1°  mean  gravity  data,  the  r.m.s.  variation  of  the  gravity  anom¬ 
alies  is  roughly  reduced  to  half  the  original  value,  and  geoid  undulations  may 
be  computed  with  an  accuracy  around  1  m  using  such  spherical  harmonic  expansions 
(Tscherning,  1983).  Thus  by  using  long-wavelength  information  from  such  expansions 
and  combining  with  short  and  medium  wavelength  topographic  effects  computed  using 
a  residual  terrain  model  with  respect  to  a  corresponding  spherical  harmonic  expan¬ 
sion  of  the  topography  itself,  the  "remaining"  signal  will  be  smooth,  its  variance 
low  and  its  degree  of  anisotropy  usually  less.  This  will  be  demonstrated  later 
in  this  report,  through  investigations  of  local  empirical  covariance  functions 
and  power  spectra  of  the  gravity  field  in  areas  with  different  types  of  topography. 

The  computation  of  terrain  effects,  using  digital  terrain  models,  is  basically 
a  problem  of  numerical  integration.  However,  it  is  by  no  means  a  simple  problem. 
The  integration  kernels  are  usually  singular  at  the  evaluation  points,  and  the 
influence  of  the  "inner  zone"  -  the  topography  in  the  immediate  vicinity  of  the 
evaluation  point-  is  very  critical  for  quantities  like  gravity  anomalies  and  second 
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order  derivatives.  A  FORTRAN  program  for  computing  terrain  effects  on  geoid  undulations 
deflections  of  the  vertical  and  gravity  anomalies,  based  on  the  rectangular  prism 
as  integration  element,  will  be  given  in  the  appendix. 

For  gravity  anomalies  a  type  of  topographic  effect  -  the  conventional  gravi¬ 
metric  terrain  correction  -  is  of  special  interest.  The  terrain  correction  is 
basical ly  not  a  terrain  reduction,  used  in  conjunction  with  a  general  gravity 
field  modelling  procedure  it  represents  no  unique  density  anomaly  distribution 
to  be  removed  from  the  observations.  Rather  the  terrain  correction  should  be 
viewed  as  a  mathematical  convenience,  representing  the  -  usually  small  -nonlinear 
part  :f  the  total  terrain  effect.  Unfortunately  terrain  corrections  have  from 
time  to  time  been  confused  with  terrain  reductions  proper.  The  application  of 
terrain  corrections  alone  does  usually  not  improve  results  of  the  gravity  field 
modelling  in  mountainous  areas  significantly,  since  the  bulk  of  the  topographic 
density  anomaly  distribution,  causing  short  wavelength  "noise"  in  the  gravity 
field,  is  not  removed.  The  application  of  terrain  corrected  free-air  anomalies 
does,  however,  make  good  sense  for  integral  formula  applications,  since  the  applica¬ 
tion  of  the  terrain  correction  to  free-air  anomalies  represents  a  first  (and  rather 
crude)  approximation  to  the  problem  of  downward  continuation  of  gravity  observations 
from  the  surface  of  the  topography  to  the  geoid,  the  terrain  correction  being 
an  approximation  to  Molodensky's  G  L  -term,  see  e.g.  Heiskanen  and  Moritz  (1967) 
and  Moritz  (1966).  The  role  of  the  terrain  correction  will  be  given  attention 
later  in  this  report,  and  some  examples  of  its  magnitude  will  be  given. 

To  summarize,  the  emphasis  in  this  report  will  be  on  the  utilization  of  density 
anomalies  in  local  gravity  field  modelling  -  especially  collocation  and  related 
methods.  The  first  part  will  review  principles  for  the  utilization  of  known  and 
unknown  density  anomalies,  then  the  practical  computation  of  such  effects  - 
especially  topographic  effects  -  will  be  outlined,  and  finally  the  influence  of 
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the  topography  will  be  studied  through  investigations  of  empirical  covariance  functions 
for  various  test  areas  in  the  United  States.  No  major  examples  of  actual  applications 
of  the  methods  for  gravity  field  modelling  will  presently  be  given.  For  earlier 
results  of  gravity  field  modeling  by  collocation  using  some  of  the  terrain  reduction 
concepts  presented,  the  reader  is  referred  to  e.g.  Forsberg  and  Tscherning  (1981), 
Forsberg  and  Madsen  (1981),  and  Tscherning  and  Forsberg  (1983). 


2.  The  Anomalous  Gravity  Field  and  Density  Anomalies 

The  gravity  field  of  the  earth  is  traditionally  described  using  the  anomalous 
potential  T 

T  =  W  -  U  (2.1) 

representing  the  difference  between  the  actual  geopotential  W  and  a  normal  poten¬ 
tial  U,  corresponding  to  chosen  reference  ellipsoid  parameters.  In  U  is  also 
included  the  centrifugal,  tidal  and  atmospheric  potentials,  and  thus  T  is  a 
harmonic  function. 


V2T  =  0 


(2.2) 


outside  the  surface  of  the  earth,  and  may  be  expanded  in  spherical  harmonics 

(2.3) 


T<r.  •■■)=?!  j 


i  -2  m=-  i 


r  P  (sin  $ )  cos  mx  (m  >  0) 
(sin  $)  sin  (m  <  0) 


Here  G  is  the  gravitational  constant,  M  the  mass  of  the  earth  and  R  the  radius 
of  a  reference  earth  sphere  (Bjerhammar  sphere). 


The  observable  gravity  field  quantities  may  in  the  usual  spherical  approx¬ 
imation  be  expressed  as  linear  functionals  L(T),  the  most  important  quantities 
being  point  and  area  mean  values  of 
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where  y  is  normal  gravity. 

Similarly,  density  anomalies  Ao  may  be  defined  as  the  difference  between 
the  actual  density  distribution  p  inside  the  earth  and  a  normal  density  distri 
bution  p  ,  generating  U: 


Ap  =  p  _  p, 


W  =J  ^dVn  +  *,  r-  |r| 
V  w 


U  E  dVQ  +  0 


(2.10 


Figure  1. 

(here  •*>  is  the  centrifugal  potential,  V  the  interior  of  the  earth  and  E  the 
reference  ellipsoid.  We  thus  have 


T(P)  =  f  Afi.  dV 
JV  r  avQ 


(2.11) 


In  other  words,  Ap  is  a  density  distribution  generating  the  anomalous  gravity 
field. 

Due  to  the  fundamental  ambiguity  of  potential  field  theory,  an  infinite  variety 
of  density  distributions  satisfying  (2.10)  exists.  If  a  spherical  normal  potential 
U  is  chosen,  indeed  any  radial  symmetric  density  distribution,  having  the  correct 
GM-value,  generates  U.  It  is  therefore  clear  that-  the  observed  gravity  field 
is  of  no  use  in  determining  a  realistic  normal  density  p  .  Instead  we  must  get 
information  on  pQ  from  other  geophysical  sources:  seismic  body-wave  travel 
times,  surface  wave  dispersion  curves,  eigen  periods  of  the  earth's  free  oscilla¬ 
tions  and  the  moment  of  inertia.  Examples  of  current  earth  models,  applicable 
for  "defining"  pQ,  is  the  HB-I  (Haddon  &  Bullen,  1969)  model  and  the  PEM-models 
(Dziewonski  et  al.,  1976). 

To  account  for  the  non-spherical  part  of  pQ  ,  we  may  resort  to  perturbing 
the  interior  density  distribution  by  small  amounts,  given  by  the  hydrostatic  equili¬ 
brium  theory.  The  flattening  of  the  interior  density  discontinuities  will  thus  be 
decreasing  downwards,  from  1/298  at  the  surface  to  1/350  at  the  core/mantle  bound¬ 
ary.  Alternatively  we  may  resort  to  a  stringent  analytic  representation  of  the 
normal  density  distributions  using  ellipsoidal  coordinates  (where  the  flattening 
increase  with  depth) ,  and  more  or  less  arbitrary  mathematical  constraints  to  secure 
a  unique  solution  (Moritz,  1968;  Tscherning &  Sunkel ,  1980).  In  any  way,  however, 
the  non-spherical  perturbations  are  very  small,  much  less  than  the  errors  in  the 
geophysical  earth  models,  and  we  may  thus  for  all  practical  purposes  simply  disre¬ 
gard  these. 


We  are  thus  free  to  choose  "convenient"  reference  density  distributions  when 


working  in  given  regions:  a  typical  continental  choice  would  be  e.g.  a  density 
starting  at  2.67  g/cm3  at  sea  level,  increasing  to  2.9  at  the  base  of  the  crust, 
jumping  to  3.3  across  the  moho  at  32  km  depth  and  increasing  through  the  mantle 
with  major  "discontinuities"  at  the  phase  transition  zones  at  ~420  km  (olivine-spinel) 
and  at  -700  km.  At  the  base  of  the  mantle  the  PEM-model  gives  a  density  of  5.4, 
and  for  the  earth's  cere  values  from  9.9  to  13.0  at  the  center,  the  density  of 
the  inner  core  being  still  very  uncertain.  For  an  oceanic  area  we  might  change 
this  model  above  the  low  velocity  zone,  e.g.  choosing  a  reference  model  with  4  km 
of  water  (density  1.0'),  a  thin,  dense  crust  (2.9)  extending  to  12-18  km  depth 
and  an  "undepleted",  cceanic  upper  mantle  at  3.4 g/cm3. 


3.  On  the  Use  of  Spherical  Harmonic  Expansions 

When  we  use  a  spherical  harmonic  expansion  as  a  first  step  in  gravity  field 
modelling,  the  "wanted"  approximation  T  to  I  is  split  into 


T  s  Ti  * T. 


with  T  given  by  the  expansion 


T,  (r,  <t> ,  \)  s  “  I 
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(3.2) 


The  currently  available  high  degree-and-order  models  (*tnax  =  180)  provides  the 
bulk  of  T.  They  suffer,  however,  of  a  minor  problem  relating  to  the  continental 
topography:  information  on  the  higher-degree  coefficients  stem  from  analysis 
of  l°x  1°  Ag  (used  directly  in  the  Rapp  models  (Rapp,  1981)  and  through  Stokes' 
formula  in  GEM10C  (Lerch  et  a!.,  1981)),  treated  as  data  on  a  sphere,  neglecting 
that  the  continental  anomalies  are  actually  anomalies  at  altitude.  This  fact 
gives  rise  to  a  small  correction,  completely  corresponding  to  Molodensky’s 
G  -term,  but  in  the  frequency  domain. 


Let  the  standard  surface  harmonic  expansion  of  the  mean  anomalies  be 
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To  first  order  these  anomalies  correspond  to  elevation  h($,  a),  defined  through 
a  similar  expansion  of  the  continental  topography  (0  at  oceans).  Using  the  correct 
spatial  expansion  (3.2)  we  get 
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where  Ag*  is  the  gravity  anomalies  harmonically  downward  continued  to  the 

Bjerhammar  sphere.  Since  a„  =  a '  ,  the  second  term  in  (3.6)  may  be  evaluated 

*  m 

with  sufficient  accuracy  from  existing  solutions,  representing  essentially 

h  •  T  .  Expanding  this  correction  term  in  surface  spherical  harmonics  b  ,  we 

tm 

obtain 
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which  by  (3.3)  and  expansion  of  Ag*  gives 

a  =  a '  +  — — b 

tm  tm  t-1  tm 


(3.7) 


(3.8) 


The  correction  term  has  for  gravity  anomalies  a  maximum  value  of  c.  19  mgal  (Rapp, 
1983).  For  local  gravity  field  modelling  the  above  has  the  practical  application 
that  elevations  of  the  individual  (ground)  observation  points  should  be  used 
properly  when  evaluating  (2.3),  otherwise  elevations  should  rather  be  set  to  zero. 


Corresponding  to  (3.1)  also  the  density  anomaly  may  be  split  in  a  spherical 
harmonic  reference  part  and  a  residual 


Ap  =  Ap l  +  Ap2 


(3.9) 


The  reference  density  distribution  Api  poses  some  problems,  especially  for  high 
degree-and-order  fields  (*  >  180),  since  many  of  the  major  crustal  -upper- 

mantle  structures  (trenches,  rifts,  etc.)  will  indeed  have  a  significant  part 
of  their  gravitational  signal  in  the  reference  part  Ap  .  When  working  with 
residuals  ("T  ")  only,  the  response  from  assumed  Ap-models  must  thus  be  split, 
either  by  introducing  a  "formal"  Ap  ,  or  by  high-pass  filtering  the  response. 

In  this  case  we  will,  however,  lose  important  information  about  the  structure. 

For  more  local  gravity  field  modelling,  we  may  totally  neglect  the  density 
split  (3.9).  Many  of  the  typical  intracrustal  density  anomalies  would  have  only 
small  long-wavelength  effects.  By  removing  such  density  anomalies  computationally. 


the  remaining  part  of  the  residual  potential  T2  would  in  principle  be  "con¬ 
taminated"  with  these  long  wavelength  errors,  but  they  will  usually  not  be  very 
significant  compared  to  e.g.  the  errors  in  the  reference  field  T  . 

For  the  topography,  the  natural  choice  of  Apt  would  be  a  model  corresponding 
to  an  analogous  expansion  in  spherical  harmonics  of  the  topographic  elevations: 


i  *max  * 

hU,  \)  =  j  |  l  V  Ytm^’  x ^ 
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(3.10) 


The  reference  density  model  in  this  case  would  have  density  -  2.67  g/cm3  below 
the  mean  elevation  surface  h"( ^> ,  \)  on  the  continents.  More  on  this  (i.e.,  the 
residual  terrain  correction)  later. 


Formal  introduction  of  spherical  harmonic  reference  density  anomalies 
may  be  done  using  simple  analytical  inversion  methods.  Consider  e.g.  a  two-layer 


2- 


earth  model  with  a  interface  at  depth  D  (Figure  2).  The  effect  of  undulations 
h (0  ,  x  )  of  this  interface  may  be  approximated  by  a  mass  coating  of  density 
<  =  ^>2  -  p^h.  We  thus  have: 


p 


Figure  2 


T(P)  =  GJ  -p^ljdofQ)  (3.11) 

1  % 

where  a  ^  is  the  interface  sphere. 


(3.11)  represents  a  spherical  convolution,  and  we  have  the  simple  well-known  ex¬ 
pression  in  the  frequency  domain  for  the  dimensionless  coefficients  (3.2)  of  the 
generated  potential 


a 


m 


4ttR2  1  /  D\* 

M  2i+l  V  ‘  R/  Km 


(3.12) 


see  e.g.  (Sunkel,  1981b).  Just  a  single  interface  thus  provides  a  unique  inversion 
For  real  applications,  however,  several  layers  will  be  needed  in  order  that  the 
derived  interface  undulations  be  reasonable  and  the  corresponding  stress  levels 
within  accepted  limits.  Typically  the  lower  harmonics  could  be  modelled  as 
"topography"  on  the  mantle/core  interface  and  the  deeper  mantle  phase  transition 
zones,  the  higher  harmonics  on  discontinuities  in  the  upper  mantle  and  the  moho. 
Such  a  model  makes  reasonable  physical  sense,  it  has  e.g.  been  hypothesized  that 
major  global  features  of  the  geoid  corresponds  to  thermally  induced  shifts  in 
the  olivine-spinel  transition  zone. 
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Alternatively,  unique  "spatial"  density  anomaly  models  may  be  obtained  by 
imposing  "analytic"  constraints  on  the  possible  density  distributions.  If  e.g. 
the  density  distribution  fulfills  the  condition 


v2(r"r)  »  0 


(3.13) 


where  n  is  an  arbitrary  integer  constant,  the  density  solution  is  found  by  an 


expansion  in  internal  spherical  harmonics  as 


(r,  ,,  X)  =  r-n[  Ibtm  (£)*  (♦.  A) 


(3.14) 


k  =  (2*  -  n  +  3)(2t  +  1)  a 
Dfcm  4^ GR^“"  urn 


(3.15) 


for  2i  >  n  -  3  (Tscherning,  1974).  The  drawback  of  this  method  is  that  the  con¬ 
dition  (3.13)  is  completely  arbitrary  without  any  physical  meaning.  The  resultant 
density  distribution  will  have  its  extremes  at  the  surface  of  the  sphere,  and  the 
actual  density  variations  will  be  very  low  -  e.g.  order-of-magnitude  0.004  g/cm3 

for  GEM10B  («■  =36)  (Tscherning  &  Sunkel,  1980).  Attempts  to  find  other  constraints 

max 

like  (3.13)  corresponding  to  some  simple  physical  minimum  principle  have  been 

fruitless  (Tscherning,  personal  communication)  -  it  is  obviously  not  possible 
to  find  "state  equations"  for  the  earth's  interior  relating  only  to  the  density 
distribution. 


4.  Utilization  of  Known  Density  Anomalies 


In  external  gravity  field  modelling  known  (or  assumed)  density  anomalies 
may  be  taken  into  account  by  a  simple  remove-restore  technique:  the  influence 
of  the  anomalous  masses  is  subtracted  from  the  given  data  ("observations"),  then 


the  gravity  field  modelling  tech¬ 
nique  is  applied  on  these  terrain 
reduced  data,  and  the  final  results 
("predictions")  are  obtained  by 
adding  back  the  terrain  effects 
to  the  predicted  anomalies. 

Let  V  be  the  volume  enclosing 
the  given  density  anomalies. 

Then  in  a  point  P 


TJP)  =  G  J  ^dVp,  r=|7g-7p|  (4.1) 

is  the  terrain  effect  potential,  and  for  a  gravity  field  quantity  L(T)  we  have 
the  "terrain"  effect  (including  "geologic"  effects) 

L(Tm)  =  G  {  Ap  1(7)  dV  (4.2) 

V  r 

e.g.  for  the  gravity  disturbance  vector 

6^m  =  _vTm  =  _G  ^ v  dV  =  G  JyA p  -fir  dV  =  [  G  Aoi  /  -jfr  dV  (4.3) 

1  vi 

For  practical  computations  "building  blocks"  of  constant  density  are  traditionally 
used,  as  expressed  by  the  last  term  of  (4.3). 

The  remove -re store  technique  may  schematically  be  written  as: 
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OBSERVATIONS: 

terrain  I 
reduction  | 


model  ling  | 
technique  | 


"inverse" 

terrain 

reduction 

PREDICTIONS: 


L.obs(T) 

Li0bs(Tc)  =  L.0bs(T)  - 


L  pred(tC) 
J 


L  pred(T) 

J 


l-jPred(fC)  *  Lj(Tm) 


Irrespectively  of  the  gravity  field  modelling  technique  actually  used  (Integral 
formulas.  Collocation  etc.),  Tc  must  be  a  harmonc  function.  This  is  secured  if 
Tm  represents  the  gravitational  effects  of  a  given,  fixed  mass  model,  e.g.  the 
density  anomalies  within  a  given  geographical  area.  The  same  mass  model  must 
naturally  be  used  for  both  the  observed  and  predicted  quantities.  For  topographic/ 
isostatic  effects  and  -  especially  -  "residual"  topographic  effects,  a  global 
mass  model  is  often  appropriate:  formally  we  thus  have  to  extend  the  integral 
(4.2)  all  around  the  earth,  but  in  practice  it  is  sufficient  to  integrate  out 
to  a  certain  distance  from  the  computation  point  -  depending  on  the  type  of  gravity 
field  quantity  and  thus  the  "sharpness"  of  the  integral  kernel  L(-p)  -  the  effect 
of  the  distant  topography  being  either  negligible  or  obtainable  from  e.g.  spherical 
harmonic  expansions. 

When  using  a  remove- res tore  technique  like  outlined  here,  it  is  important 
to  know  that  the  assumed  density  anomalies  need  not  be  realistic  -  an^  density 
distribution  may  be  used  as  long  as  Tm  and  thus  Tc  is  harmonic  outside  the  topography. 
But  naturally  the  most  smooth  Tc  is  expected  when  the  most  realistic  mass  model 
is  applied.  For  "geologic"  density  anomalies  -  e.g.  salt  domes,  intrusions,  faults 
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etc.  -  even  the  simplest  models  (spheres,  cylinders,  step  functions)  may  often 
be  applied  successfully,  to  give  a  more  stationary  and  isotropic  residual  field, 
well-suited  for  interpolation  and  approximation. 

5.  Unknown  Densities  -  Geophysical  Inversion 

Most  frequently  we  do  not  have  a  good  knowledge  of  "geologic"  density 
anomalies,  and  it  would  therefore  be  natural  to  try  to  estimate  parameters  des¬ 
cribing  such  density  anomalies  -  preferably  simultaneously  with  the  external 
gravity  field  modelling  process.  In  addition  to  the  obvious  importance  of  know¬ 
ledge  of  the  surface  density  distribution,  we  will  by  this  method  also  have  a 
simple  way  of  introducing  non-gravity  observation  data  (magnetic,  seismic  etc.) 
in  the  external  gravity  field  modelling.  Simple  but  very  essential  geodetic 
applications  includes  determination  of  optimum  topographic  reduction  densities 
(generalized  Nettleton  method)  and  e.g.  determination  of  ocean  bathymetry  in  un¬ 
surveyed  areas  from  satellite  altimetry. 

The  general  principles  of  geophysical  inversion  nay  be  outlined  as  follows: 
based  on  geologic  intuition  (or  practical  ease),  a  (finite)  number  of  parameters 
xi ,  i  =  l,...n  is  chosen  to  represent  the  structures  in  a  given  area  (Fig.  4). 

A  gravity  field  quantity  observed  in  a  point  P  may  then  formally  be  expressed 


Figure  4 


as 


Lp(T)  =  fp{xi . xn) 


(5.1) 


•k  - 


(5.2) 


where  the  function  fp  is  generally  non-linear  and  must  be  linearized 

(Lp(T)  -  Lp(T) )  =  l  ^  ax. 

i  Q 

by  assuming  an  initial  model  x?.  Since  geophysical  inversion  problems  are  often 
highly  nonlinear,  a  large  number  of  iterations  (5.2)  are  often  necessary.  The 
model  parameters  x^  may  be  generally  classified  in  two  types: 

1)  geometric  parameters  (interface  depths  etc.) 

2)  density  value  parameters 

The  main  emphasis  in  traditional  geophysical  modelling  has  been  in  terms 
of  structural  geometric  parameters  (see  e.g.  Burkhard  &  Jackson,  1976;  Pedersen, 
1979),  to  directly  represent  interfaces  such  as  the  top  basement  in  sedimentary 
basins  or  the  moho,  exemplified  in  Figure  4  (left).  The  advantage  of  the  geo¬ 
metric  parameters  is  that  they  directly  reflect  simplified  geologic  models,  and 
additional  data  such  as  well  control  is  easily  including  by  e.g.  fixing  one  or 
more  parameters.  The  drawback  of  choosing  "structural"  parameters  is  the  inherent 
unlinearities. 

Opposed  to  this,  models  with  density  value  parameters  only  (Figure  4,  right) 
are  perfectly  linear,  but  the  computational  advantage  of  the  linearity  is  usually 
counterweighted  by  the  greater  number  of  parameters  needed  to  represent  a  wanted 
geologic  scenario.  Also  it  is  less  simple  to  include  the  non-gravity  constraints. 
The  commonly  applied  point  mass  modelling  in  geodesy  may  be  viewed  as  a  special 
case  of  such  geophysical  inversion,  using  the  simplest  possible  finite  element 
representation  (delta  spikes)  of  the  subsurface  density  distribution.  However, 
this  simplest  possible  case  of  inversion  gives  results  that  are  amazingly  close 
to  results  obtained  using  improved  (spatial)  density  representations  (Figure  5). 


Assuming  a  set  of  observations  y.  =  L.(T),  i  =  l,...,n  the  (linearized) 
inversion  geophysical  problem  may  be  written  as 


y  =  A< 


(5.3) 


where  A  is  the  response  matrix.  This  problem  is  generally  ill-conditioned  or 
improperly  posed,  and  generalized  invers  on  must  not  be  used. 

One  popular  technique  is  the  use  of  the  singular  value  decomposition: 


(5.4) 


with  the  U  and  V 
ATAv  =  x2v 

J  J  J 

AATu  «  x?u. 

J  J  J 


being  orthonormal  matrices  defined  through 
V  =  tvj} 

U  =  (uj( 


(5.5) 


See  e.g.  Pedersen  (1979)  or  Rummel  et  al.  (1979).  p  is  the  number  of  non-zero 
eigen  values,  i.e.  the  number  of  degrees  of  freedom  of  the  problem.  A  solution 

x  to  (5.3)  is  given  by  the  Lanczos  inverse, 

x  =  VA'Vy  (5.6) 

minimizing  as  well  y^y  as  x^x.  To  prevent  the  eigen  values  of  ill-conditioned 
problems  to  induce  large  changes  in  the  parameters  x,  the  eigen  value  spectrum 

xi  may  be  truncated  by  removing  eigen  values  smaller  than  a  suitable  threshold, 
giving  the  traditional  trade-off  between  resolution  and  variance. 

Alternatively  to  the  explicit  use  of  the  singular  value  decomposition,  es¬ 


sentially  the  same  solution  may  be  obtained  using  Tikhonov  regularization.  In 
this  case  we  seek  to  minimize  a  combination. 


Assuming  a  noise  covariance  matrix  D  for  the  observations  and  an  a  priori  covariance 
matrix  C  for  the  "signal"  x  (both  matrices  usually  assumed  to  be  diagonal),  we 
obtain  the  solution  by  solving  the  normal  equations. 

(ATD_1A  +  aC"1 )x  =  ATD_1y  (5.8) 

• 

(Runnel  et  al . ,  1979).  The  constant  a  is  arbitrary  and  may  be  chosen  to  obtain 
a  desired  smoothness  of  the  parameters,  again  with  the  price  to  be  paid  being 
a  degraded  fit  of  the  model. 

Independent  geologic  information  may  be  taken  into  account  using  linear  con¬ 
straints  of  the  form 

Bx  =  C  (5.9) 

where  B  and  C  are  constant.  Such  constraints  can  be  used  to  fix  certain 
parameters  (e.g.  representing  known  depths  to  an  interface),  to  fix  differences 
in  density  values  (e.g.  forcing  parameters  of  type  "2"  to  represent  layers,  faults, 
etc.)  and  to  introduce  special  geometric  constraints  on  the  anomalous  mass  body 
based  on  geologic  experience  (e.g.  assuming  a  dike  to  have  parallel  sides). 

The  constraint  (5.9)  is  taken  into  account  in  the  minimization  problem  (5.8)  using 
Lagrange  multipliers,  obtaining  somewhat  more  complicated  normal  equations. 

Details  may  be  found  e.g.  in  Burkhard  &  Jackson  (1976). 

The  methods  outlined  above  represent  conventional  geophysical  inversion  tech¬ 
niques.  They  are  usually  applied  only  for  one  type  of  gravity  field  quantity 
(gravity  anomalies  or  -  at  times  -  altimeter  geoid  undulations),  but  there  is 
of  course  no  restriction  in  the  model  formulation  to  utilize  heterogeneous  data 
(e.g.  simultaneous  gravity  and  geoid  information)  as  we  are  commonly  used  to  in 
geodesy.  The  problem  with  the  heterogeneous  data  lies  in  the  reqional/residual 
separation:  .the  gravity  field  cont. ins  information  about  density  anomalies  at 
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all  depths,  but  the  model  parameters  x  are  typically  restricted  to  describe 
simplified  rather  shallow  structures  -  a  filtering  is  therefore  done  to  remove 
the  unwanted  parts  of  the  signal.  This  filtering  is  often  very  crude  (e.g.  graph¬ 
ical  determination  of  a  "regional")  and  not  applicable  for  heterogeneous  data, 
for  such  data  we  must  make  sure  that  the  filtering  of  the  different  data  types 
are  consistent  -  the  "regional"  must  be  a  harmonic  function. 

In  some  cases  high  degree  and  order  spherical  harmonic  expansions  might  be 
valuable  as  “regionals"  -  e.g.  when  trying  tc  model  total  crustal  density  distri- 
butions-but  we  should  then  also  have  a  well-defined  spherical  harmonic  reference 
density  distribution  (c.f.  Section  3).  Alternatively  we  can  utilize  "general" 
gravity  field  modelling  techniques  to  represent  the  regional,  e.g.  by  introducing 
arbitrary  (deep)  model  point-masses  or  by  doing  the  inversion  within  the  framework 
of  least  squares  collocation  with  parameters. 

In  this  case  we  have  the  following  observation  equations  for  an  observation 
yi  with  associated  linear  functional  ard  noise  n^: 

yi  *  0\xl.  +  L.(T)  +  n.  (5.10) 

for  which  we  get  the  well-known  collocation  solution  (see  e.g.  Moritz,  1980) 

T(Q)  -  LjK( • ,Q)  T  C"1 

x  =  (ATC“lA)“:  ATC-1 y  (5.11) 

C  =UiU(v)^D.j) 

where  D  again  is  the  noise  covariance  matrix,  K(P,  Q)  the  potential  covariance 
function  of  the  gravity  field.  Note  that  this  covariance  function  should  not 
be  the  observed,  empirical  covariance  function  but  rather  the  covariance  function 
of  the  field  after  the  model  influence  have  been  subtracted  -  i.e.  the  covariance 
function  of  the  "regional".  We  would  expect  this  field  to  have  less  variance 
and  greater  correlation  length  than  the  original  field.  Since  the  model  results 


depend  on  the  covariance  parameters,  these  must  ultimately  be  determined  through 
trials  or  through  considerations  of  the  wanted  characteristics  of  the  regional/ 
residual  filter. 

Least  squares  collocation  with  parameters  will  be  especially  well-suited 
for  the  determination  of  optimum  topographic  reduction  densities  in  mountainous 
areas.  In  this  case  our  model  parameters  x  will  just  be  a  single  value  (or 
a  few,  if  the  geology  is  changing),  and  the  observation  equation  (5.10)  will  look 
1  i  ke 

yi  =  {G  ^  Li<7)dV}  Ap  +  MT)  +  ni  (5.12) 

where  the  term  in  the  bracket  represents  the  terrain  effect  of  a  topography  with 
unit  density,  cf.  (4.2).  This  prob  em  is  well-conditioned  for  sufficiently 
varying  topography,  and  represents  \  straight  forward  generalization  of 
Nettletons  density  profiling  method  to  heterogeneous  data.  More  reliable  density 
estimates  are  obtained  with  (5.12)  :han  with  the  more  traditional  approaches  such 
as  regression  analysis  of  the  varia:ion  of  free-air  anomalies  with  elevation, 
as  pointed  out  by  Sunkel  (1981a).  Application  of  (5.12)  will  probably  be  even 
better  than  using  real  measirements  of  sample  rock  densities:  everybody  who  has 
tried  this  knows  how  difficilt  it  is  to  estimate  average  formation  densities  from 
samples  of  individual  rock  formations,  especially  for  sedimentary  rocks  with  their 
varying  porosity  and  water  saturation. 

When  estimating  more  complex  structural  models  of  the  density  anomalies, 
stabilization  of  the  parameters  x  in  (5.11)  will  be  needed,  and  we  will  have 
to  make  a  combined  collocation/jeneral ized  inverse  approach.  Collocation  by  itself 
may  be  viewed  as  an  inversion  problem  (Moritz,  1976):  the  simple  collocation 
approximation  T  is  built  up  from  the  kernel  function  K(P,  Q)  in  the  observation 
points: 
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T(Q)  =  l  ajLjK(  •  >Q)  (5.13) 

J 

where  the  coefficient  a.  is  the  solution  fo  the  "normal  equations"  corresponding 

j 

to  (5.11).  Expressing  (5.6)  in  terms  of  these  coefficients  we  have: 


■ *x  + 1  lil/  ■  aj + ni 


-  (a  !  {117} +  "1 


(5.14) 


which  clearly  shows  our  problem  as  a  "double"  generalized  inverse  problem  with 
unknowns  x.  .  (geophysical  parameters)  and  a.  (kernel  coefficients).  The  solution 

J  J 

is  obtained  by  minimizing  a  combination: 

II  x  I!2  +a||T||2  (5.15) 

where  a  is  a  positive  constant  and  ||  ||H  the  Hilbert  space  norm  associated  with 
the  chosen  covariance  function  K.  The  constant  a  is  arbitrary,  and  must  be 
chosen  based  on  empirical  investigations.  The  constant  determines  how  much  vari¬ 
ation  is  put  "into"  the  structure  and  how  much  is  retained  in  the  outer,  residual 
field,  and  acts  like  the  "trade-off"  parameter  in  (5.7).  By  combining  the  well-known 
methods  of  collocation  and  geophysical  generalized  inversion  like  outlined  here, 
we  have  in  fact  obtained  a  discrete  version  of  the  so-called  "mixed  collocation", 
suggested  by  Sanso  and  Tscherning  (1982). 

The  practical  applicability  of  hybrid  gravity  field  modell ing/geophysical 
inversion  methods  remains  to  be  seen.  For  geodesy  and  external  gravity  field 
modelling  the  obvious  application  would  lie  in  the  determination  of  only  a  few 
key  parameters:  topographic  densities,  density  contrasts  across  major  known  discon¬ 
tinuities  (e.g.  for  moho  at  continental  margins)  and  density  anomalies  of  well-known 
geologic  bodies  (e.g.  salt  domes),  avoiding  unlinear  structural  parameters  requiring 
iteration.  The  computational  burden  would  not  be  significc ntly  increased  using 
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such  a  limited  set  of  parameters,  and  by  choosing  "good"  geologically  reasonable 
parameters,  one  could  hope  in  many  cases  to  get  significant  improvements  in  the 
characteristics  of  the  "background"  field:  less  power,  more  stationarity  and 
a  higher  degree  of  isotropy. 

Probably  the  geophysical  exploration  would  benefit  more  from  the  hybrid  col¬ 
location/inversion  scheme.  With  the  technological  advances  heterogeneous  gravity 
field  data  will  be  more  common  -  through  the  development  of  high-precision  inertial 
survey  systems  measuring  the  complete  gravity  vector,  through  airborne  gradiometry 
and  through  geoid  undulations  from  GPS  and  satellite  altimetry  in  addition  to 
terrestrial  or  airborne  gravity.  To  perform  quantitative  interpretations  with 
error  analysis  etc.  for  such  data,  some  kind  of  "hybrid"  inversion  method  will 
be  necessary,  to  stringently  handle  model  oversimplifications,  regional/residual 
separations  etc. 

With  these  remarks  the  general  discussion  of  density  anomalies  and  inversion 
techniques  will  be  concluded.  In  the  next  section  formulas  for  actual  computations 
will  be  given, and  then  the  main  density  anomaly  -  the  topography  -  will  be  treated 
in  detail . 

6.  Density  Modelling  Using  Rectangular  Prisms 

6.1  Space  Domain 

For  the  practical  evaluation  of  gravitational  effects  of  density  anomalies, 
integrals  of  the  type: 

L(T)  =  G  J  ApL(-^)dV  (6.1) 

m  V  ” 

must  be  computed  numerically.  This  computation  is  most  naturally  done  using  the 
simplest  form  of  finite  element  representation  of  the  density  distribution:  assuming 
the  density  anomaly  Ap  to  be  constant  in  subblocks,  each  such  finite  element 


(subblock)  being  a  rectangular  prism.  For  terrain  reductions  using  digital  models, 
these  subblocks  e.g.  naturally  correspond  to  the  subdivision  defined  by  the  eleva¬ 
tion  data  grid.  The  evaluation  of  integrals  (6.1)  over  each  finite  element  is 
synonymous  with  the  formulas  for  the  gravitational  effects  of  the  rectangular 
prism  of  constant  density. 


To  integrate  spherica  symmetric  function  like  over  an  interval  with 
Cartesian  symmetry  is  doomed  to  give  some  very  complicated  integrals,  this  being 
indeed  the  case  for  the  rectangular  prism  formulas.  Let  the  coordinate  system 
used  have  axes  parallel  to  the  prism  sides  and  origin  in  the  computation  point 
P,  as  indicated  in  Figure  (>.  In  the  sequel  r  =  (x,  y,  z)  is  the  coordinate  of 
the  integration  point  Q  n  this  system.  We  have  in  P  for  various  gravimetric 
quantities: 


<o  y,  z. 


GApf  -JdV  =  GApf  /  / 
v  r  Zi 

~  dxdydz,  r  =  /“x2+yz+z2 

(6.2) 

-GH  3l^)dV  ’  'H 

— ~  (-jr)dV  =  -Gap/  fdV 
oZq  r  y  r 

(6.3) 

- -,:P  WQ 

(4)  .  -GAo|  dV 

vr3/  y  rs 

(6.4) 

■  <*»>  ■  "*>  V 

dV 

(6.5) 

Since  differentiations  occur  under  the  integrals  for  the  higher  order  derivatives, 
these  will  give  the  simple. t  formulas.  Let  the  formulas  (6.2)  -  (6.5)  be  evaluated 


as  undefinite  integrals  to  keep  the  notation  simple.  We  have  then  for  the  second 
order  derivatives 


Tzz:  -Gap  /  \  jy  dxdy  =  -GApz  /  £dy  =  GAp  z  arctan  (|£) 

x  y  y 

T  :  GAp  /  /  -K  dxdy  =  GAp  J  -£■  dy  =  GAp  log  (y+r) 
yz  r  y r 


(6.6) 

(6.7) 


For  the  first  order  derivative  a  non-trivial  integration  of  (6.7)  with  respect 
to  x  is  obtained  (Jung,  1961): 


6 g :  -Gap  J  /  y  dxdy  =  Gap  f  log  (y+r)dx  = 

—  xy  x 

GAp ( x  log  (y+r)  +  y  log  (x+r)  -  z  arctan  (|^))  (6.8) 

Finally  the  formula  for  geoid  undulations  (height  anomalies)  are  obtained  by  inte¬ 
grating  (6.8)  with  respect  to  z  (MacMillan,  1958): 


T:  GAp[xy  log  (z+r)  +  xz  log  (y+r)  +  yz  log  (x+r) 

-  y  arctan  ^  ^  arctan  yy  -  y  arctan  |^] 


(6.9) 


The  final  formulas  for  the  rectangular  prisms  are  obtained  by  summing  the 

expressions  (6.6)  -  (6.9)  over  the  corners  of  the  prisms  with  alternating  signs, 

2  2  2  . 

e.g.  T  =  l  l  l  (-l)1+J+kT..u.  where  T  .  is  (6.9)  evaluated  at  (x. ,  y,,  z.). 

i=i  j=i  k=i  1J  1JK  i  J  k 

The  formulas  for  the  remaining  derivatives  (deflections  of  the  vertical,  other 
second  order  gradients)  are  simply  obtained  by  coordinate  permutations,  see  Forsberg 
and  Tscherning  (1981). 

Although  some  simplifications  of  the  final  formulas  are  possible  using  addition 
theorems  for  logarithms  and  arctan  (arctan  a  +  arctan  b  =  arctan  a.  -  hb  ) ,  the  formulas 
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are  still  very  complex  and  time  consuming.  In  the  terrain  effect  computation 
program  (see  appendix)  an  approximative  formula,  where  the  mass  of  the  prism  is 
condensed  as  a  mass  layer  on  the  xy-plane  through  the  center  of  the  prism,  is 
used  for  geoid  undulations  instead  of  (6.9).  In  this  case  we  get  an  integral 
similar  to  (6.8): 


T  =  Gap (z2  -  z , )  \  f  ^  dxdy 

X  y  z=zm 
m 

=  Gap(z2-z1)  !|  x  log  (y+r)  +  y  log  (x+r) 


2m  =  £iiLL, 


(6.10) 


z  arctan  — 
m  z  r 

m 


yz 

yi 


For  terrain  effect  computations,  this  formula  has  negligible  error  (typically 
corresponding  to  millimeters  in  c). 

In  larger  distances  from  the  prism,  the  formulas  (6.6)  -  (6.9)  may  be  sub¬ 
stituted  by  much  simpler  series  expressions  of  the  gravity  field,  obtained  using 
a  spherical  harmonic  expansion  of  the  prism  field.  Since  the  spherical  harmonics 
expressed  in  cartesian  coordinates  are  simple  homogeneous  polynomials  in  x,  y,  and 
z,  the  resulting  series  expansions  are  simple.  In  a  prism-centered  coordinate 
system  we  have  for  the  potential 

T  =  Gap  AxAyAzj-p  +  -^rr  [(2Axz-Ayz-Az2)xz  +  (-Axz+2Ayz-Azz)yz 

\  (6.11) 

+  (-AX2-Ayz+2AZ2)zz]  +  28gr9  [oCjX4  +  a2y4  +  . . .]  +  ...  j  , 

AX  =  X2-Xl,  Ay=y2-y,,  az=z2-Zi 

(MacMillan,  1958),  from  which  gravity,  second  order  derivatives  etc.  are  easily 
found  by  differentation.  The  first  term  in  (6.11)  is  simply  the  point  mass 
approximation.  The  second  term  takes  into  account  the  different  dimensions  of 
the  prism  -  it  is  zero  for  a  cube. 


*•  -  “  •  ■ 
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In  the  terrain  effect  computation  programme  given  in  the  appendix,  a  sub¬ 
routine  "PRISM"  forms  the  nucleus  of  the  calculations.  This  subroutine  uses  the 
exact  prism  formulas  when  the  computation  point  P  is  near  the  prism,  in  an  inter¬ 
mediate  zone  the  MacMillan  formula  (6.11)  is  used,  and  finally  in  very  large  dis¬ 
tances  the  point  mass  approximation  is  used.  The  shift  between  the  formulas  is 
automatic,  determined  by  accuracy  levels  wanted  by  the  user.  (cf.  Figure  7). 

It  is  through  the  additional  use  of  the  approximate  formulas  that  the  prism  method 
becomes  feasible  for  "routine"  gravity  field  modelling  in  mountainous  areas,  other¬ 
wise  evaluation  of  the  complex  exact  prism  formulas  would  often  become  prohibitive 
in  terms  of  computer  time.  Furthermore,  in  large  distances  the  formulas  (6.6)-(6.10) 
become  numerically  unstable,  requiring  extended  precision  due  to  the  differencing 
of  the  evaluated  "corner"  -  functions  entering  the  formulas. 

6.2  Frequency  Domain 

While  the  prism  formulas  are  complicated  in  the  siace  domain,  they  are 
surprisingly  simple  in  the  frequency  domain.  Since  th»  basics  of  Fourier  analysis 
of  potential  fields  is  not  generally  well  known  in  geoiesy,  a  short 
outline  will  be  given  first. 

The  Fourier  spectral  analysis  is  applied  in  the  flat-earth  approximation. 

Let  tt  be  the  reference  plane  (e.g.  sea  level)  with  c)ordinates  (x,  y),  and 
the  associated  spectral  plane  with  spatial  frequencies  (u,  v).  Then  the  Fourier 
transformation  is  given  by  denotes  transformed  quintities): 

T(u,  v)  =  {  T(x,  y)e_i ^ux+vy^dxdy  (6.12) 

7T 

T(x,  y)  =^7  /  T(u,  vje1  ^ux-*’vy^dudv  (6.13,' 

TT 

Upward  continuation  of  the  field  to  elevation  z  is  obtained  by  a  filtering 


T(u,  v,  z)  =  T(u,  v)e"wZ 


U)  =/u^+V2 


(6.14) 
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Figure  7  Approximate  maximal  approximation  error  between  the  prism  formulas  ((6.8), 
(6.10))  and  the  simpler  MacMillan  and  point  mass  formulas  (6.11).  Errors  given 
in  percent  as  a  function  of  normalized  distance  to  prism  (r)  and  height  (h)  for  a 
square  sector,  i.e.  graphs  show  errors  in  computed  terrain  effects  from  a  rec¬ 
tangular  mountain,  with  unit  side  lengths  and  height  h  in  distance  r  from  the 
center.  For  a  cube  (h=l)  the  MacMillan  and  point  mass  formulas  are  identical.  Note 
that  for  the  geoid  the  comparison  is  against  the  "mass  plane"  formula  (6.10).  The 
graphs  are  intended  as  a  guide  for  deciding  the  accuracy  of  the  terrain  effect 
computation  program  (apoendix),  which  essentially  uses  the  value  of  r  to  discrim¬ 
inate  between  the  various  formulas.  (Example:  if  topography  is  given  on  a  1000  m 
grid  with  elevatons  up  to  2000  m  (h=2),  a  maximal  1%  error  requires  r-5,  i.e. 
the  MacMillan  formula  can  be  used  for  topography  more  than  5  km  away  from  the 
computation  point). 


and  similarly  for  the  gravity  field  functionals  (2.5)-(2.7)  simple  linear  filters 
transform  the  quantities  in  the  spectral  domain: 


!(u,  v) 

=  -4fT(u,v) 

(6.15) 

ft ( U ,  v) 

-  t  v) 

(6.16) 

A~g(u,  v) 

2  - 

=  (-W  -  -p)  T(u,  v)  =  -u)T(u,  v) 

(6.17) 

where  R  in  (6.17)  is  the  radius  of  the  earth. 

For  radial  symmetric  functions,  f(x,  y)  =  f(r'),  r '  =  / x2+yz ,  the  Fourier 
transform  (6.12)  becomes  a  Hankel  transform  (Papoulis,  1968): 

f(u,  v)  =  2?rf (id )  to  =  /u2+vz  (6.18) 

where  the  Hankel  transform  pair  (transform  denoted  by  a  bar)  is  given  by: 

f(u)  =  f°  r'  f(r')  J  Ur' )  dr'  (6.19) 

o  o 

f(r')  =  J00  u  f(«)  Jn  Ur')  dw  (6.20) 

o  u 

Here  JQ(-)  is  the  Bessel  function  of  order  zero.  Of  special  importance  is  the 
Hankel  transform  of  the  inverse  distance: 


1  _  1  Hankel  i  _,„7 

r'TF^P  « - - 


(6.21) 


(Papoulis,  1968,  p.  145). 

Now,  for  a  rectangular  prism  (Figure  6),  situated  below  the  x-y  plane,  we  have 
T(x,  y,  0)  =  Gap  [  ^  dx'dy'dz'  , 


(6.22) 


giving  the  transform 


T(u,  v)  =  Gap  /  j i.  e_1  (ux+vy)  dx'dy 'dz 'dxdy 

71  V 

■  2ttGao  f  —  e"“z  e~^ux'+vy '^dx’dy'dz' 

‘\l  U 

=  2ttGap  (e  '°Z2-e  JZl)  e~i(ux+vy)  *2  y?-  (6.23) 

Results  for  gravity  and  deflections  may  be  obtained  from  (6.23)  using  (6. 15 ) - ( 6 . 17) 
and  using  (6.18)  and  (6.21)  by  interchanging  the  order  of  integration.  Formulas 
like  (6.23)  have  been  used  for  a  number  of  years  in  geophysical  exploration,  especially 
for  the  magnetic  field  (Bhattacharyya ,  1964). 

The  advantage  of  the  formula  (6.23)  is  that  it  allows  the  use  of  the  fast 
fourier  transform  (FFT)  when  computing  the  gravity  effects  from  a  regular  grid  of 
prisms,  e.g.  defined  through  a  digital  terrain  model.  If  we  have  a  set  of  nxm 
prisms,  the  corners  of  the  prisms,  projected  on  the  x-y  plane,  will  form  a 
(n+l)x(m+l)  grid  mesh.  By  rearranging  the  sum  (6.23)  as  sums  over  this  grid,  the 
general  expression  for  the  total  effect  of  all  prisms  will  have  the  form 

f(u,  v)  =  n+l  T  f(»,  x.,  y  )e"i(uXJ'+vyk)  (6.24) 

j=l  k=l  J  K 

where  f  contains  sums  and  differences  of  e"™2  for  the  prisms  adjoining  the  grid 
point  (x.,  y  ).  Sums  like  (6.24)  is  exactly  what  is  obtained  by  the  FFT  algorithm  - 

J  K 

had  it  not  been  for  the  dependence  of  f  with  w.  This  dependence  is  due  to  the 
basic  fact  that  the  prism  integral  (6.22)  is  fundamentally  uni  inear,  not  being  a 
convolution.  We  are  thus  forced  to  evaluate  (6.24)  on  a  frequency-by-frequency 
basis  by  FFT,  for  each  value  of  w  a  separate  spectrum  T  is  obtained  and  the 
final  spectrum  must  then  be  "constructed"  by  carefuly  selection  and  interpolation 
in  this  set  of  spectra.  The  thus  obtained  final  spectrum  may  then  be  transformed 
back  into  the  space  domain  by  an  inverse  FFT. 


It  is  important  to  stress,  that  (6.24)  is  exact.  Therefore  the  spectral 
values  obtained  using  (6.24)  are  not  influenced  by  window  effects  etc.,  the  ob¬ 
tained  spectrum  represents  the  spectrum  of  a  transient  signal,  this  'signal  decreasing 
quickly  to  zero  outside  the  area  covered  by  the  prisms.  The  only  "errors"  occurring 
in  this  FFT  technique  is  in  the  ^-interpolation  scheme  to  obtain  the  final  spectrum, 
and  in  the  final  synthesis  of  the  frequencies,  since  FFT  only  gives  the  sums  (6.24) 
for  a  finite,  discrete  number  of  frequencies,  the  highest  frequencies  being  the 
Nyquist  frequencies  for  the  prism  grid.  This  secures,  however,  a  nice  smoothness 
of  the  computed  field,  since  e.g.  a  representation  of  the  topography  with  flat- 
topped  prisms  is  anyway  a  rather  poor  representation,  causing  unwanted  high  frequency 
spectral  "ripple"  effects  from  the  edges. 

To  estimate  the  gain  in  computing  speed,  consider  as  an  example  a  nxm  grid 
of  prisms  (with  varying  top  and  base  levels),  and  assume  we  want  to  compute  the 
gravitational  effects  in  the  same  grid  at  a  fixed  altitude.  Then  the  operations 
will  be  (orders  of  magnitude): 

SPACE  DOMAIN:  n2xn2  calls  of  "PRISM"  subroutine  (no  computation¬ 

saving  grid  symmetries  exists  for  exact"  formulas) 

FREQUENCY  DOMAIN:  n/2  spectra  (6.24)  of  n2  coefficients  f, 

FFT  speed  -  n2logN,  spectral  selection,  inverse  FFT. 

Combined  order  of  magnitude:  n3logN 

The  gain  is  thus  moderate,  a  consequence  of  the  unlinearity  of  (6.22). 

A  real  significant  gain  in  computation  speed  is  obtained  if  the  basic  volume 
integral  (6.1)  is  approximated  with  surface  convolution  integrals.  Inis  is 
e.g  possible  for  "thin"  prism  layers  at  near  constant  depth,  and  to  some  degree 
also  for  terrain  effects  (so-called  "linear  topographic  approximation"),  involving 
integrals  of  the  topographic  elevations  and  their  squares  (more  details  in  next 
section).  In  the  case  of  a  "thin  layer"  at  average  depth  D,  we  obtain 


(6.25) 


T( x,  y,  0) 


G  J  dx'dy'dz ' 
V  r 


G  f 


*  (x1 ,  y1 ) 


"  [(x-x'^+ty-y'^+D2]^. 


-  dx'dy' 


where  <  =  Ap  ( z2-z1 )  is  the  surface  density.  The  transform  is  obtained  simply 
by  utilizing  (6.21)  again,  giving 


T(u,  v)  =  2^G-  e-a)D  k  (u,  v) 

Ui 


(6.26) 


In  this  case  the  order-of-magnitude  computation  speed  of  the  previous  example 
will  be  only  n2logn  if  FFT  is  utilized,  but  opposed  to  the  "exact"  spectral  formu¬ 
lation  "window  effects"  due  to  finite  data  lengths  must. now  be  given  full  attention. 

The  frequency  domain  methods  have  as  common  restrictions  that  data  and  compu¬ 
tation  points  must  be  in  a  grid,  the  computation  points  being  in  a  plane  (im¬ 
portant  exception:  grtvimetric  terrain  corrections,  cf.  next  section).  Obvious 
applications  could  be  e.g.  for  geoid  computations  at  sea  level  (especially  for 
satellite  altimetry)  and  upward  continuation  studies  (airborne  gravimetry  and 
gradiometry) .  The  importance  of  spectral  methods  in  geophysical  inversion  may 
be  inferred  from  (6.26):  if  a  particular  spectrum  (e.g.  white  noise)  is  expected 
for  the  source  <,  then  the  depth  D  to  the  sour:e  may  be  found  directly  from 
the  observed  gravity  field  spectrum.  This  is  the  base  of  the  widespread  "statis¬ 
tical  inversion  techniques",  dominating  in  the  analysis  of  aeromagnetic  data. 


7.  Terrain  Reductions 

For  the  remainder  of  this  report,  emphasis  will  now  concentrate  on  topo¬ 
graphic  and  isostatic  reductions  -  a  synonym  for  computational  elimination  of 
the  effects  of  the  two  most  dominant  and  best  known  density  anomalies  of  the  earth: 
the  visible  topography  and  its  associated  compensation  at  depth.  For  such  gravity 
field  effects  the  general  term  "terrain  effects"  will  be  used  in  the  present  context 
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The  commonly  applied  term  "terrain  corrections"  will  lie  reserved  for  the  narrow 
meaning,  i.e.  a  correction  to  the  Bouguer  reduction,  to  give  the  true  (unlinear) 
effect  of  the  topography  on  gravity  anomalies  (and  deflections  of  the  vertical 
as  well). 


7.1  Terrain  Effects  and  Associated  Density  Anomalies 


The  various  terrain  reductions  in  use  is  illustr.ted  in  Figure  8.  To  use 
terrain  reduced  data  in  a  "remove-restore"  technique  ior  gravity  field  modelling 
(Section  4),  it  should  be  remembered  that  the  density  models  indicated  by  Figure 
8  should  either  cover  a  given,  fixed  geographical  area,  or  -  or  at  least  in  prin¬ 
ciple  -  be  global . 


The  topographic  reduction  or  complete  Bouguer  reduction  consists  of  removing 
the  visible  topography.  Conventionally  a  density  of  2.67  g/cm3  is  used.  This 
density,  which  represents  a  typical  density  of  granite  and  many  Paleozoic  and 
Pre-Cambrian  sediments,  is  fairly  good  in  mountainous  areas.  However,  one  should 
not  hesitate  to  use  other  density  values,  since  the  density  may  range  from  below 
2.0  g/cm3  in  moraine  hills  to  3.0  g/cm3  in  volcanics.  Average  density  values 
could  be  chosen  from  geological  considerations  or  usirg  the  inversion  techniques 
of  the  last  section.  At  the  oceans  the  topographic  density  anomalies  are  formally 
negative,  the  standard  density  2.67  correspondi ng  to  1 .03-2.67= -1.64  g/cm3,  1.03 
being  the  density  of  sea  water. 

The  topographic  reduction  may  formally  be  split  into  a  Bouguer  term,  the 
effect  of  an  infinite  plate,  plus  the  terrain  correction,  which  takes  into  account 
the  topographic  irregularities.  For  gravity  disturbances  we  have 


<59topo  =  2lTGohP  '  tC  (7-15 

where  27rGphp  is  the  gravity  due  to  a  (plane)  Bouguer  plate  of  thickness  hp 
(if  the  computation  point  P  is  situated  above  the  topography,  hp  is  the  topo¬ 
graphic  elevation  at  the  surface  point  below  P),  and  tc  is  the  terrain  correction. 
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Figure  8.  Density  anomalies  associated  with  various  terrain  reductions  (con¬ 
tinental  area).  A:  topographic  effect,  i.e.  the  "complete"  Bouguer  reduction 
(consisting  of  the  effects  of  a  Bouguer  plate  minus  the  terrain  correction  "B"), 
C:  conventional  Airy-isostatic  model,  D:  Residual  terrain  model  (RTM),  the 
mean  elevation  surface  e.g.  given  by  a  180  x  130  spherical  harmonic  expansion. 
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The  terrain  correction  is  always  positive  (in  the  plane  approximation)  due  to 
the  conventional  minus  sign  in  (7.1). 

For  deflections  of  the  vertical  the  terrain  correction  and  the  topographic 
effect  are  identical  (except  the  signs),  since  the  Bouguer  plate  effect  is  zero. 

For  height  anomalies,  however,  the  infinite  Bouguer  plane  can  not  be  used,  the 
effect  being  infinite.  Instead  one  could  think  of  using  a  spherical  Bouguer  plate: 
the  effects  typically  computed,  will  however,  still  be  very  large  and  often  much 
larger  than  the  observed  geoid  undulations  themselves  (on  a  global  basis).  Topo¬ 
graphic  reductions  are  therefore  not  very  applicable  to  general  gravity  field 
modelling:  the  large  model  geoid  effects  and  biased  Bouguer  anomalies  at  oceans 
and  mountainous  areas  necessitates  some  kind  of  negative  density  anomalies  being 
introduced,  e.g  through  an  isostatic  compensation  hypothesis.  Needless  to  say, 
the  topographic  reduction  is  naturally  very  well  suited  for  problems  such  as  gravity 
interpolation  and  geophysical  inversion. 

Isostatic  reduction  formalizes  the  prevailing  tendency  of  the  earths  topo¬ 
graphy  to  be  compensated  at  depth.  The  standard  Airy  scheme  assumes  local  compen¬ 
sation  through  a  root  system  (Figure  8D),  the  thickness  of  the  root  being 

t  *  ^  h  !  6.7  h  (7.2) 

where  o  is  the  density  of  the  topography  (~2.67  g/cn3)  and  Ao  the  density 
contrast  between  the  crust  and  the  mantle  (-0.4  g/cm3).  The  normal  density  model 
has  a  crust  of  thickness  D  (~  32  km). 

Naturally  the  earth  does  not  fully  follow  this  simple  principle.  Although 
(7.2)  approximates  the  overal 1  isostatic  compensation  fairly  good,  many  exceptions 
occur:  first  of  all  the  strength  of  the  earth's  crust  supports  short-wavelength 
topographic  features,  isostasy  being  primarily  a  regional  phenomena.  Second, 
many  regions  are  either  uncompensated  or  compensated  at  deeper  levels  (through 


anomalous  density  values  in  the  upper  mantle) , most  noticeably  the  deep-sea  trenches 
and  mid-oceanic  ridges.  However,  since  the  computed  isostatic  effects  are  very 
insensitive  to  the  actual  isostatic  formulation  and  parameters  used,  even  the 
simplest  formulations  (e.g.  (7-2))  gives  excellent  results,  the  results  being  "good" 
when  the  remaining  field  after  isostatic  reduction  is  smooth  and  with  low  variance. 
Global  isostatic  reductions  attain  maximal  values  for  the  geoid  in  the  range  10-20  m 
It  is  therefore  necessary  to  compute  isostatic  effects  also  on  spherical  harmonic 
coefficients  .for  the  geopotential,  e.g.  using  the  simple  formula  (3.12). 

A  drawback  of  the  isostatic  reduction  is  that  it  primarily  should  be  global. 

If  only  a  fixed,  localized  area  is  taken  into  account,  the  computed  isostatic 
gravity  field  effects  will  be  influenced  by  "edge  effects":  the  computed  isostatic 
gravity  and  deflections  of  vertical  would  change  rapidly  near  the  boundary  for 
non-zero  elevations.  For  the  geoid  an  overall  bias,  dependent  on  the  chosen  size 
of  the  reduction  area,  will  result  if  the  area  mean  elevation  is  different  from 
zero  (see  e.g.  Forsberg  &  'scheming,  1981,  Figure  1). 

Since  the  main  problem  in  external  gravity  field  modelling  in  mountainous 
areas  is  short-wavelength  topographic  "gravity  field  noise",  a  terrain  reduction 
method  avoiding  the  "global"  computations  of  isostatic  reductions,  but  capable 
of  approximating  isostatic  conditions,  would  be  ideal: 

For  a  residual  terrain  model  (RTM)  reduction  only  the  short  wavelength 
of  the  topography  is  taken  into  account.  This  is  done  by  choosing  a  smooth  mean 
elevation  surface,  and  computationally  remove  masses  above  this  surface  and  fill 
up  valleys  below  (Figure  8D).  The  mean  elevation  surface  could  be  any  smooth 
surface,  representing  mean  elevations  of  the  area,  e.g.  an  interpolation  in  30 ' x30 ' 
mean  heights  or  -  especially  -  defined  through  a  high-order  spherical  harmonic 
expansion  of  the  topography  of  the  earth.  In  this  case  the  RTM  density  anomalies 
correspond  to  a  normal  density  distribution  (normal  earth)  with  smooth  topography 
and  bathymetry  defined  through  the  spherical  harmonic  expansion,  and  thus 
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corresponds  to  the  residual  gravity  field  after  removal  of  a  similar  spherical 
harmonic  expansion  of  the  geopotential. 

The  advantages  of  the  RTM-reduction  are  many:  snce  the  density  anomalies 
have  oscillating  positive  and  negative  values,  the  iniegrations  for  gravity  field 
effects  need  only  be  done  out  to  some  suitable  distance  from  the  computation  point, 
the  influence  of  distant  topography  cancelling  out.  /Iso,  terrain  effects  on 
height  anomalies  will  be  small  (often  negligible  if  a  short-wavelength  reference 
elevation  surface  is  chosen),  and  especially  for  e.g.  180  x  180  height  expansions 
the  reduction  gives  results  close  to  isostatic  reductions. 

The  similarity  between  RTM  and  isostatic  reductions  are  analogous  to  the 
similarity  between  mean  free-air  gravity  anomalies  ana  isostatic  anomalies. 

Indeed,  by  a  special  choice  of  mean  elevation  surface  nearly  complete  correspon¬ 
dence  may  be  obtained:  If  we  define  the  mean  elevations  through  the  low-pass 
filter  (plane  approximation) 


:(P)  / 


2-  ;  [r2+D2]3/2  • 


r  =  dist  (P,  Q) 


(7.3) 


then  Moritz  has  shown  that  the  associated  RTM-reduction  corresponds  to  an 


isostatic  reduction  with  a  (surface  layer)  compensation  depth  D  (Moritz,  1968a). 
Note  that  (7.3)  is  nothing  but  the  well-known  Poisson  integral  for  upward  con¬ 
tinuation  of  harmonic  functions. 

The  RTM-reduction  may  be  viewed  as  a  difference  between  two  Bouguer  reductions: 
first  the  visible  topography  is  removed,  and  then  the  smoothed  topography  is  added 
back  (Figure  9): 


l(t)rtm  =  l(tHopo  "  L^T)REF-T0P0 


(7.4) 


Each  term  in  (7.4)  may  formally  be  split  in  a  Bouguer  plate  effect  and  a  terrain 
correction.  Table  1  shows  sample  terrain  corrections  for  a  180  x  180  spherical 
harmonic  reference  surface  in  two  4°x  4°  fixed  areas  in  tne  Rocky  Mountains. 
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Table  1.  R.M.S.  and  absolute  maximal  terrain  corrections  for  a  180  x  180  spherical 
harmonic  reference  topography  (4°x4°  fixed  area,  9  sample  points). 

From  the  table  it  is  seen  that  the  gravity  reference  terrain  corrections  are 
very  small  -  below  1  mgal*.  We  may  therefore  for  gravity  anomalies  simply  state 

AgRTM180  =  27rGp  (h”href^  '  tc  ^7*5^ 

i.e.,  when  using  a  RTM  reduction  with  180  x  180  reference  heights  (RTM180)  the 
terrain  effect  is  simply  a  terrain  corrected  (tc)  Bouguer  reduction  to  the  level 
hre^.  This  has  the  important  practical  advantage  that  available,  terrain-corrected 
Bouguer  anomalies  (being  still  the  bulk  of  the  available  local  gravity  field  data) 


♦Additional  verification  on  actual  data:  see  Section  7.4. 
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may  be  applied  directly  for  RTM-reduction  using  (7.5).  For  deflections  of 
the  vertical,  however,  the  time-consuming  "prism"-integrations  can  not  be  avoided, 
the  deflection  terrain  effects  due  to  the  changing  reference  level  being  much 
too  large. 

When  performing  the  RTM-reduction  "directly"  (e.g.  using  rectangular  prism 
integration),  stations  above  the  reference  level  are  left  "hanging  in  the  air", 
while  observations  below  this  level  are  reduced  to  their  values  inside  the  mass 
(Figure  9).  However,  for  external  gravity  field  modelling,  we  need  not  the  value 
inside  the  mass,  but  the  harmonically  downward  continued  value,  corresponding 
to  the  outer,  "reduced"  field.  In  other  words,  what  would  the  reduced  observation 
be  at  the  point  P2  in  Figure  9,  if  we  treated  the  mean  topography  as  non-existent 

An  approximate  answer  to  this  question  is  simple:  if  the  density  above 
a  plane  through  P2  is  condensed  in  a  mass  plane  layer  immediately  below  P2  , 
deflections  of  the  vertical  and  geoid  undulations  would  remain  nearly  unchanged 
due  to  the  smooth,  low-slope  reference  surface.  For  gravity  anomalies,  however, 
we  would  see  a  change 

Agharmonic  '  A^in  mass  ~  4irGpAh  (7.6) 

corresponding  to  a  "double"  Bouguer  reduction  with  plate  thickness  Ah  =  h^  -  hp. 
This  "harmonic  correction"  must  be  applied  for  all  gravity  stations  below  the 
reference  level  when  "direct"  prism  integration  of  RTM  density  anomalies  (Figure 
8D)  is  performed.  If  instead  (7.5)  is  used,  the  correction  is  taken  into  account 
"implicitly" . 

7.2  Practical  Terrain  Reductions  in  Gravity  Field  Modelling 

A  FORTRAN  77  program  for  computation  of  any  of  tha  four  types  of  terrain 
effects  (and  corrections)  mentioned  (Figure  8)  are  listed  in  the  appendix. 


The  program  uses  rectangular  prisms  for  a  direct  integration  of  geoid  undulations, 
deflections  of  the  vertical  or  gravity  anomalies  from  digital  terrain  models  given 
on  a  geographic  grid. 

Special  precautions  have  been  taken  to  evaluate  the  inner  zone  effect,  i.e. 
the  influence  of  the  topography  adjacent  (say,  within  1  km)  to  the  computation 
point.  These  inner  zone  effects  may  be  very  large,  especially  for  gravity  terrain 
corrections.  To  represent  the  inner  zone,  a  bicubic  spline  interpolation  of  the 
topography  is  utilized.  However,  since  gravity  topographic  effects  to  first  order 
depends  linearly  on  the  gravity  station  elevation,  it  is  clear  that  the  station 
elevation  itself  should  be  utilized  in  the  inner  zone  interpolation.  An  option 
in  the  program  allows  the  height  interpolation  procedure  to  give  the  correct  eleva¬ 
tion  at  a  station .through  a  smooth  "adjustment"  of  the  digitial  terrain  model 
elevations  in  the  inner  zone.  For  deflections  and  height  anomalies,  where  the 
station  height  dependence  is  weak,  use  of  this  option  is  not  necessary. 

Actual  examples  of  use  of  the  various  terrain  reductions  in  connection  with 
gravity  field  modelling  by  collocation  can  be  found  in  e.g.  Forsberg  and  Tscherning 
(1981).  Here  gravity  and  deflections  were  modelled  with  an  accuracy  around  4mgal 
ana  1"  respectively,  in  a  mountainous  area  (New  Mexico),  using  gravity  data  spaced 
c.  6'  apart  and  a  0 . 5 ' jc 0.5*  digital  terrain  model.  When  applied  properly,  nearly 
the  same  results  were  obtained  for  all  types  of  terrain  reductions. 

As  an  outline  example,  let  us  consider  upward  continuation  of  gravity  data 
in  a  mountainous  area.  Using  a  "spatial"  modelling  technique  like  collocation 
or  point  mass  modelling  the  appl ication  of  the  remove-restore  technique  for 
a  RTM180-reduction  (and  a  180  x  180  reference  field)  is  simple: 

1.  Compute  terrain  corrections  for  local  gravity  stations  if  not  already  given 

2-  Obtain  terrain-reduced  residual  gravity  values  by  subtracting  the 

reference  Bouguer  anomalies  Ag^£p  -  2^GohRgp  from  the  local,  terrain- 
corrected  Bouguer  anomalies. 

3.  Apply  upward  continuation  method, 

4.  Add  back  RTM-effects  computed  at  altitude  (prism  integration), 

5.  Add  back  180  x  180  gravity  computed  at  altitude. 
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For  already  gridded  gravity  data  (e.g.  5 ‘ x  5 '  mean  free-air  anomalies)  this 
remove- restore  technique  may  be  used  with  some  caution.  For  a  mean  block  we  would 
need  the  mean  terrain  correction,  as  we  have  from  (7.5) 

AgRTM180  =  21x2(3  ^  “  ^REF^  '  tc 

Such  mean  terrain  corrections  tc  are  difficult  to  estimate.  They  are, 
however,  very  important  since  they  play  an  essential  role  in  the  harmonic  down¬ 
ward  continuation  of  gravity  data  from  the  surface  of  the  topography  to  the  geoid, 
a  necessary  prerequisite  for  the  application  of  e.g.  the  classical  integral  methods. 
Apart  from  direct  computation  of  tc  by  averaging,  its  magnitude  may  be  estimated 
from  the  covariance  function  of  the  topography 

-  -  °h 

tc  =  3tt  Gp  -jj1- 

where  a ^  is  the  terrain  variance  and  d  the  correlation  length,  as  pointed 
out  by  Sunkel  (1981a) . 

7.3  The  Linear  Approximation  for  Topographic  Effects 

Approximate  formulas  for  RTM-effects,  especially  applicable  for  error  studies 
and  frequency  domain  methods,  may  be  obtained  using  functional  expansions  of  the 
topographic  volume  integral  kernels  (y  ,  yj  etc.).  In  the  sequel  a  "long  wave¬ 
length"  reference  elevation  surface  e.g.  180  x  180  spherical  harmonic  expansion 
is  assumed  to  be  used. 

In  the  plane  approximation  we  have  for  the  RTM  potential  effect  when  a  con¬ 
stant  topographic  density  is  used: 

h  . 

Tm  =  G  j  ^-dV  =  Go  J  J  j  dzd*  (7.8) 

*  z=href 


where 
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7  =  tr0  +  (z-*)2rH  *  7-[l  =  -i-[l  +  tan2e Th  (7.9) 

r  0  p  ro  r  02  ro 

V  is  the  volume  of  non-zero  (positive  or  negative)  density  anomalies  Ao »  shown 
hatched  in  Figure  10.  A  series  expansion  for  Tm  is  obtained  by  expanding  (7.9) 
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Figure  10 


with  respect  to  the  inclination  8,  which  is  always  small  except  for  the  inner 
zone  in  rugged  topography  (and  geoid  innerzone  effects  are  very  small).  Thus, 

1  =  J_  .  1  (2  "  hp)2  +  ...  (7.10) 

r  r  2  r  3 
o  o 

and  by  inserting  (7.10)  in  (7.8)  and  integrating  with  respect  to  z 
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(7.11) 


The  first  term  in  (7.11)  has  been  called  the  linear  approximation  by  Moritz  (1968a) 
It  is  seen  that  this  term  represents  nothing  but  the  potential  of  a  mass  coating 
<  =  p(h  -  href)  =  pAh.  Within  the  accuracy  of  the  linear  approximation  we  may 
view  this  mass  coating  as  a  surface  density  layer  at  the  height  reference  surface. 
The  higher  order  terms  in  (7.11)  may  similarly  be  viewed  as  successive  coatings 
of  dipoles,  quadropoles,  etc. 


For  deflections  of  the  vertical  we  have 


cm  r  y  -  yD  .  r  h  cos  a 

„  i-i  f  p  =  j  / 

m  j  v^x-Xpl  1  it  z=hrep  I  sin  a 


(7.12) 


which 


by  expansion  r-3  =  rg3  -  ^  rjj5  (z  -  hp)2  +  ...  and  integration  analogous 


to  the  potential  case  gives 


,Gp_  f  Lh^refj  (C0Sa)H^lGL  f  (h-hp)3  -  (href-hp)3 
y  1  r2  2  Y  * 


(7.13) 


Aqain,  the  first  term  (the  linear  approximation)  may  be  interpreted  as  the  effect 
of  a  mass  coating. 


The  RTM  gravity  anomalies  are  given  by 


A9m  =  -Gp  J  /  ^7^  Tm 


TT  Z  =  href 


(7.14) 


The  last  term  in  (7.14)  -  the  indirect  effect  -  will  usually  be  below  1  mgal 
for  a  180  x  180  reference  surface  and  may  be  neglec  ed.  For  gravity  anomalies 
it  is  advantageous  to  use  (7.5) 


Agm  =  2uGp  (h  -  href )  -  tc 


and  then  only  expand  the  terrain  correction 


(7.15) 


tc  =  Gp  j '  | 


z  -  h , 


(7.16) 
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In  (7.15)  the  Bouguer  term  is  on  the  average  one  order  of  magnitude  larger  than 
the  terrain  corrections  themselves.  In  some  cases  -e.g.  error  studies  -  it  would 
therefore  be  sufficient  only  to  use  this  simple  term.  For  gravity  field  modelling 
with  heterogeneous  date  it  is,  however,  very  dangerous  not  to  include  the  best 
possible  terrain  corrections:  since  tc  is  always  positive,  a  systematic  bias 
will  be  introduced  in  egm,  a  bias  which  often  would  seriously  affect  computed 
geoid  undulations. 

For  frequency  domain  formulas  we  note,  that  the  linear  approximations  in 
(7.11)  and  (7.13)  are  convolutions,  i.e.  expressions  of  the  form 

L(Tm)p  =  k*Ah  =  J  k(xp-Xg,  yp-yg)  Ah  (xQ,  yg)dirQ  (7.18) 

and  thus  in  principle  for  the  Fourier  transform  using  the  well-known  convolution 
theorem 

L(Tm)  -  k  Ah  (7.19) 

The  function  k  is  traditionally  called  the  impulse  response,  k  the  transfer 
function.  Ah  =  h-href  will  be  termed  the  residual  height  hereafter. 

For  the  potential  we  have  from  (6.26) 

Tm(u,  v)  =  2ttGo  Ah  (u,  v)  ,  <*>  *  /  u?+vz  (7.20) 

1  sT 

and  since  (also)  in  the  linear  approximation  f  -  -  —  ®  we  have 

m  Y  ojr 

and  by  (7.15) 


Ag  (u,  v)  =  2nGp  Ah  (u,  v)  -  tc  (u,  v) 


(7.22) 


Note,  that  (7.20)  and  (7.21)  may  be  derived  from  the  Bouguer  term  of  (7.22), 
excluding  the  terrain  correction,  however,  we  might  < xpect  the  linear  approximation 
for  potential  and  deflections  to  be  "better"  than  the  Bouguer  aproximaticn  for 
gravity:  intuitively  the  condensation  of  the  "rod"  oi  Figure  10  to  a  point  at 
the  mean  elevation  surface  would  have  a  large  effect  cn  gravity  but  nearly  none 
on  potential  and  deflections.  More  formally,  from  the  expansions  of  the  impulse 
responses  in  terms  of  the  inclination  6  : 

potential  and  deflections:  +  k383  +  ... 

gravity:  ’  kQ  +  k2B2  +  k4B4  +  ... 

it  is  seen  that  the  "condensation"  interpretation  (Boiguer  term  for  gravity)  cor¬ 
responds  to  a  first  order  expansion  in  8.  By  includ'ng  the  terrain  correction 
for  gravity  anomalies,  a  second  order  expansion  in  8  is  obtained.  It  is  thi s 
expansion  which  Moritz  (1968a)  has  termed  the  "linear"  approximation ,  in  order 
to  have  the  same  accuracy  level  as  in  the  well-known  linearized  Molodensky  approach 
to  the  geodetic  boundary  value  problem. 

7,4  Accuracy  of  the  Linear  Approximation 

If  the  linear  approximation  is  sufficiently  accurate,  much  faster  techniques 
for  the  evaluation  of  terrain  effects  in  grids  are  available  (FFT  and  space  domain 
filtering  methods). 

As  a  simple  analytical  example,  consider  the  terrain  correction  at  the  summit 


of  a  cone  shaped  mountain. 


The  exact  terrain  correction  at  P  is: 


h  _  Li  i 

tc  =  Gp  /  /  dzd tt  =  Gp  /  ( 

7T  H  r 


IT  o 


-)  dTT 


(7.23) 
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2tt gpH  sine 


(7.24) 


while  the  linear  approximation  gives 


2  ~r  J  r- 

IT  0 


diT 


=  1  Go  f  (rQ-  9-)  -  2rrrQdro  +  j  Gp  /  -p-  2irrndr 
0  0  so 


0  0 


2irGpH  tane  (7.25) 


The  relative  error  for  some  slope  values  e  are: 


0 

tc 

Agm  (total  topographic  effect) 

15° 

3.5% 

1.2% 

15% 

15% 

45° 

41% 

100% 

For  common  slopes  the  linear  approximation  thus  seems  somewhat  reasonable,  since 
the  cone  mountain  is  a  "worst-case"  model. 


For  a  practical  evaluation  of  the  linear  approximation  (and  the  error  associated 
with  the  "Bouguer-split"  (7.5)  for  gravity  anomalies),  an  alpine  l°x  1°  block  in 
Colorado  has  been  chosen  (Figure  12).  Comparisons  have  been  performed  in  36 
stations,  located  in  a  12 * x  12'  grid  at  the  surface  of  the  topography.  As  eleva¬ 
tion  data  0. 5 ' x  0. 5'  heights  were  used ,  covering  totally  a  4°x  4°  area  surrounding 
the  comparison  area.  To  get  the  linear  approximation  results,  the  "prism"  sub- 
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routine  of  the  terrain  effect  computation  program  (cf.  appendix)  was  changed  to 
give  instead  only  mass  pi ane  results  (for  formulas  see  Forsberg  &  Tscherning, 
1981).  For  the  linear  approximation  to  the  terrain  correction,  we  have  for 

a  prism  element  at  position  Xj  to  x2,  yx  to  y2  relative  to  the  computation  point: 


tc  -  -jjGp 


*2  y 2 

f  f 
X1  *1 


(h-hg)2 


dxdy 


r  =/  x2+y2 


(7.26) 


which  by  simple  integration  gives 


tc  =  |  Gp  (h-hp)2 


V  1  £  .  .y* 

/  x7  r  dx 


»  -  \  Gp  (h-hp)2 


r 

2 

xy 

X, 

The  following  results  were  obtained: 


(7.27) 


Table  3.  Comparison  "exact"  RTM-reduction  versus  linear  approximation, 
Colorado/Mt.  Evans  area: 


Quant i ty 

Exact 

Mean 

Std . 
Dev. 

Linear  Approximation 
Std. 

Mean  Dev. 

Mean 

Difference 

Std. 

Abs. 

Max. 

Cm(meter) 

0.82 

0.65 

0.82 

0.65 

0.00 

0.01 

0.04 

Cm(arcsec) 

0.68 

6.85 

0.60 

6.59 

0.09 

0.69 

-1.47 

nm(arcsec) 

-0.21 

7.06 

-0.34 

6.78 

0.13 

0.72 

2.79 

Agm(mgal ) 

-7.45 

43.61 

-7.74 

43.49 

0.29 

0.60 

3.24 

tc  (mgal ) 

6.01 

_ 

4.41 

6.33 

4.94 

-0.29 

0.60 

-3.24 

Agm  in  Table  3  are  computed  using  (7.5).  A  comparison  between  a  rigorous  RTM 
prism  computation  of  Agm  and  (7.5)  gave  as  a  result  an  r.m.s.  difference  of 
only  0.3  mgal  (maximal  value  0.7  mgal),  thus  supporting  the  simple  "Bouguer"  ■ 
formula  (7.5). 


From  Table  3  it  is  seen  that  the  linear  approximation  errors  are  insignificant 
for  geoid  undulations  and  also  rather  small  on  the  average  for  gravity,  but  with 
a  possibility  for  rather  large  outliers.  For  deflections  the  error  might  not 
always  be  acceptable,  and  higher  order  expansions  might  be  necessary.  Considering 
the  extreme  ruggedness  of  the  test  area,  other  "milder"  areas  will  be  expected 
to  give  better  results,  and  certainly  the  linear  approximation  will  always  be 
very  useful  since  it  allows  the  use  of  FFT  methods  for  terrain  effect  computations. 

7.5  The  Terrain  Correction  as  Convolution  Integrals 
From  (7.17)  we  have  in  the  linear  approximation 

tc  =  hGp  !  V2-  dir  (7.28) 

*  ro 

which  may  be  expressed  as  convolutions  of  h  and  h2: 

tc  =  JjGp [  [  — —  dir  +  f  — P-  dir  -  f  ;P-  dir]  (7.29) 

*  r|  ’  r3o  -  ri 

since  hp  is  constant  with  respect  to  the  integration,  we  obtain 

tc  =  hGp  [(h2*f )  +  h2  /  f  d*  -  2h  (h*f)]  (7.30) 

”  TT  P 

where  f  =  — — — - -~sr-  .  Now,  the  function  f  does  not  have  a  fourier  spectrum, 

(x2+y2)/2 

but  it  may  be  regularized  very  simply:  Instead  of  f  consider  f'  =  - - - - jr- 

(xz+y2+a2)  /2 

where  a  is  a  small  constant.  Using  f'  instead  of  f  as  kernel  in  (7.28),  this 
corresponds  to  a  computational  upward  continuation  to  a  distance  a,  and  when 
a  is  chosen  sufficiently  small  the  error  will  be  insignificant  and  only  affect 
an  innerzone  roughly  of  radius  a.  We  thus  have: 


tc  =  J^Gp  [(h2*f 1  i  +  ^h2  -  2h  (h*f ' )] 

a  p  p 


(7.31) 


where  the  center  integral  of  (7.30)  has  been  evaluated  analytically  (the  integral 
is  nothing  but  the  well-known  Bouguer-plate  integral). 
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The  two  convolutions  in  (7.31)  may  with  advantage  be  evaluated  using  fre¬ 
quency  domain  methods,  we  have  e.g. 

transform  of  (h2*f‘)  =  h2  f  =  e"'*>a  h2  (7.32) 

see  Papoulis  (1968),  p.  145.  With  the  split  of  the  terrain  correction  (7.31) 

it  is  clear,  that  for  practical  applications  there  will  be  a  lower  limit  for  "a", 

* 

since  the  terrain  correction  is  expressed  as  a  (small)  difference  between  two 
large  numbers,  and  thus  unstable  numerically. 

Since  "a"  can  not  be  chosen  arbitrarily  small,  a  small  error  Ate  is  made, 
by  using  integral  kernel  f1  instead  of  f,  rt presenting  the  "suppressed"  effect 
of  the  local  innerzone  just  around  the  computation  point.  If  the  regularization 
distance  "a"  is  chosen  somewhat  smaller  than  the  finest  resolution  of  the  given 
elevation  data,  a  quantitative  estimate  of  tiis  "regularization  error"  Ate  may 
be  obtained  as  follows:  The  error  is 


(7.33) 


Since  this  error  is  dominated  by  a  very  local  innerzone  effect,  a  Taylor  expansion 

o 

of  the  topography  is  adequate,  keeping  only  the  first  terms  to  represent  a  sloping 
plane.  Then  by  assuming  the  computation  point  to  be  at  origo  we  have 


(h-hpr  =  r02  tan2e  cos2a 


(7.34) 


where  e  is  the  slope  of  the  plane  and  a  the  azimuth  from  the  direction  of 
maximal  slope.  Insertion  of  (7.34)  into  (7.33)  gives 

/  t*  (£irri  •  dro 


=  %ttGo  tan2e  {  (1  - 

n 


[rz+az]? 


tJ  d 


r 

o 


(7.35) 
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This  integral  may  be  tackled  by  substitution  of  r^  ,  giving 

2  2  i  oo 

Ate  =  JjirGp  tan2e  rn  -  rQ+  = .  I 

“  'o 


=  irGp  tan2ea  (7.36) 

This  is  exactly  the  (linear  approximation)  terrair  correction  of  a  sloping  innerzone 
of  radius  2a  (it  is  evident  from  (7.34)  that  this  terrain  correction  must  be  half 
the  corresponding  cone  terrain  correction  (7.25)).  A  numerical  example:  if  a  =  100  m 
and  e=  30°,  the  regularization  error  will  be  1.8  mgal.  This  illustrates  the 
critical  importance  of  the  very  local  station  surroundings  for  gravity  terrain 
corrections,  and  the  fact  that  "a"  can  not  be  chosen  too  large  -  reasonable  values 
must  be  chosen  based  on  empirical  investigations. 

7.6  The  Use  of  FFT  for  Terrain  Effect  Computations 

In  the  previous  sections  frequency-domain  formulas  for  terrain  effects  on 
height  anomalies,  deflections  of  the  vertical  and  gravity  anomalies  have  been 
discussed.  These  expressions  are  very  useful  since  digital  terrain  models  are 
naturally  given  in  grids,  suitable  for  direct  use  of  the  Fast  Fourier  Transform 
(FFT).  The  speed  of  FFT  (~nlog2n  if  the  number  of  points  n  is  a  power  of  2, 
otherwise  somewhat  slower  depending  on  the  prime  factorization  of  n)  certainly 
makes  the  application  of  frequency  domain  methods  attractive  especially  when  large 
volumes  of  data  meed  to  be  terrain  reduced.  However,  before  applying  the 
method,  it  is  essential  to  realize  the  limitations  inherent  when  using  FFT. 

FFT  is  basically  a  fast  algorithm  to  determine  the  discrete  transform  of 
periodic  data.  The  two-dimensional  discrete  Fourier  transform  pair  may  be  expressed  a 

n  —  1  m~ 1  .  ,Pj  ,  qk  » 

h(pAu,  qAv)  =  — ^  l  l  h  (j'ax,  kAy)e"/n  'n  m  (7.37) 

nm  j=0  k=0 


-  I/-  m 

.  1^  . 


■  * 
i  ■ 


W." 
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,  n-1  m-1  . 

h( jAx,  kAy)  =  j——  £  l  h(pAU,  qAv)e2  (7.38) 

p=0  q=0 

where  ax  and  Ay  are  gridspacing  in  the  given  elevation  grid  of  nxm  points, 
and  the  normalization  factors  have  been  chosen  to  be  in  accordance  with  the  con¬ 
tinuous  transform  (6.12-6.13).  The  frequency  spacings  are  given  by 


AU 


2tt 

nAx  ’ 


AV 


2~n 

mAy 


(7.39) 


The  spectrum  h  as  well  as  the  original  data  h  must  be  viewed  as  infinitely 
periodically  extended  in  space.  Since  the  Nyquist  frequencies 


n  it  m  .  it 

JN  =  2  Au  =  ax  ’  VN  "  2  Av  '  Ay 


(7.40) 


represent  the  highest  frequencies  obtainable  from  the  gridded  data,  frequencies 
above  (u^,  v^)  in  (7.37)  will  correspond  to  negative  frequencies.  For  more  details 
on  FFT  see  e.g.  Kanasewich  (1975). 

When  applying  FFT  for  convolutions  of  the  form  (7.18)  or  (7.31),  the  periodic 
extension  means  that  for  a  data  point  near  an  edge,  the  convolution  will  actually 
"use"  data  from  around  the  opposite  edge  as  well.  The  convolution  kernel  will 
in  effect  be  truncated  (and  periodically  extended)  when  applying  the  analytical 
terrain  effect  filters  at  the  discrete  frequencies  of  the  FFT  (Figure  13),  and 
this  means  in  common  words  that  the  terrain  effects  computed  using  the  FFT  method 
will  be  terrain  effects  from  a  "running  area"  of  size  (nAx,  mAy),  centered  at  the 
computation  point.  At  the  center  point  of  the  grid  the  computed  effect  will 
exactly  be  the  effect  of  the  given  area  -  at  a  corner  point  the  result  will  be 
completely  erroneous,  since  3  of  the  4  quadrants  around  the  corner  will  be  integrated 
with  the  "non-existent"  periodically  extended  heights. 


1 


elevations 


integral  kernel 
Figure  13 


terrain  effects 


For  terrain  reductions  in  connection  with  general  gravity  field  modelling, 
it  has  earlier  been  stressed  that  to  preserve  consistency  and  harmonicity  either  a 
fixed  mass  model  must  in  principle  be  taken  into  account,  or  -  expecially  for 
residual  terrain  reductions  -  the  terrain  effect  integrations  must  be  carried 
out  to  a  sufficient  distance  from  the  computation  points,  so  that  the  influence 
of  the  remote  zones  will  be  negligible  for  all  quantitites  (gravity,  deflections, 
height  anomalies  etc.).  The  FFT  methods  will  in  principle  only  be  applicable 
for  general  terrain  reductions  when  either  sufficiently  large  areas  of  elevation 
data  surrounding  the  target  area  is  transformed,  or  the  given  elevation  grid  is 
extended  with  a  "border"  of  zero-values  on  all  sides  (of  "width"  ^nax  and  %m\y 
respectively)  to  obtain  a  "true"  fixed  area  reduction  at  the  price  of  a  quadrupling 
of  the  elevation  data.  In  both  cases  computer  limitations  in  storage  might  be 
prohibitive.  Consider  e.g.  RTM-reduction,  with  a  180x180  spherical  harmonic 
reference  surface.  If  (as  is  very  common)  1  km  x  1km  elevation  data  is  available, 
an  area  of  dimension  300-400  km  must  be  taken  into  account  for  complete  reduc¬ 
tion  of  deflections  and  height  anomalies  in  a  l°x  1°  block,  necessitating  a  complex 
array  of  size  "region"  1440  k-2560  k  in  double  precision  IBM  FORTRAN.  Furthermore, 
when  the  interest  is  concentrated  in  a  rather  small  area,  like  the  l°x  1°  area, 
it  seems  somewhat  unnecessary  to  take  into  account  every  tiny  topographic  irregular¬ 
ity  at  large  distances,  which  is  in  principle  done  in  the  simple  FFT  approach. 

These  drawbacks  may  be  overcome  by  using  a  "hierarchial"  set  of  FFT  terrain 
reductions,  utilizing  more  and  more  coarse  mean  elevation  grids.  In  effect,  the 
terrain  computations  are  split  into  various  "wavelength  bands".  This  split  may 
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be  done  either  in  the  space  domain  or  the  frequency  domain.  Consider  the  simple 
case  of  two  elevation  grids,  a  detailed  DTM  (e.g  lxl  km)  covering  the  target 
area  with  a  rather  small  margin,  and  a  coarse  DTM  (e.g.  10  x  10  km  mean  heights) 
covering  a  much  larger  area. 

In  the  space  domain  the  terrain  effect  convolution  kernel  f  may  be  split  in 
two  parts,  a  near-zone  and  a  far-zone  effect,  symbolically 

terrain  effect  «•  f*Ah  =  (fj  +  f2)  ★  Ah  =  f  *  Ah  +  fo  *  Ah  *(7.41) 

fj  being  the  integral  from  zero  to  a  certain  distance  d,  f  the  integral  from 
d  to  infinity  (Figure  14).  (For  gravity  terrain  corrections  two  convolutions 
of  h  and  h2  are  needed,  as  discussed  in  the  last  section).  By  choosing  a  suitable 
d,  f2  may  be  computed  by  sufficient  accuracy  from  the  coarse  elevation  grid  and 
interpolated  to  the  points  of  the  detailed  elevation  grid,  where  the  innzer  zone 
contributions  are  evaluated  from  f1#  The  drawback  of  the  space  domain  split  is 
that  the  "truncated"  transfer  functions  fx  and  f2  do  not  have  simple  analytical 
expressions,  and  a  numerical  transform  must  be  made  to  obtain  f,  and  . 


where  the  elevation  Ah  is  split  into  a  smooth  mean  elevation  surface  Ah2»  e.g. 
obtained  by  interpolation  in  the  coarse  elevation  grid,  and  a  residual  height 
Ahf  -  analogous  to  the  RTM- reduction.  In  other  words,  the  terrain  effect  is  split  . 
into  a  long-wavelength  part  from  the  mean  elevations  and  a  short-wavelength  part 
from  the  detailed  elevations.  Ideally,  the  Ah2-surface  should  be  a  low-pass  filtered 
version  of  Ah,  so  that  Ah;  contains  only  frequencies  above  the  Nyquist  frequencies 
for  the  coarse  grid,  and  Ah2  only  frequencies  below.  Otherwise  errors  due  to 
aliasing  wi 1 1  occur.  This  may  be  illustrated  as  follows: 

Let  Ah2  be  defined  through  a  running  average: 


,  x+*sax2  y+isAy2 
Ah2 ( x ,  y)  =  —  /  /  Ah(x\  y')  dx'dy' 


^  ^ y 


.(7.43) 


where  (ax2,  Ay2)  are  the  grid  spacings  of  the  coarse  grid.  Then  Ah2  may  be 
expressed  as 


Ah; 


1 


Ax2Ay2 


g*Ah, 


g  = 


I  x|  <  Jjax2,  I  y I  <  %Ay2 

otherwise 


(7.44) 


which  can  have  energy  at  all  frequencies  since 


g(Uj  v)  =  1  sin  ( 2  u/un)  sin  (2v/vN) 


(7.45) 


where  (u^,  v^)  are  the  Nyquist  frequencies  (7.40)  for  the  coarse  grid.  (See  e.g. 
Papoulis,  1968).  Now,  if  FFT  is  directly  applied  on  the  coarse  (averaged)  heights, 
the  non-zero  spectrum  of  Ah2  above  (u^,  v^)  will  by  (7.44)  and  (7.45)  result 
in  a  non-zero  spectrum  of  Ah2  above  this  interval,  which  "folds"  erroneously 
into  the  low  frequencies  ("aliasing")  by  FFT,  to  give  long-wavelength  errors  in 
the  coarse  terrain  effect.  In  addition,  minor  errors  occur  when  interpolating 
from  points  in  the  coarse  grid  to  the  dense,  detailed  elevation  points.  The  aliasing 
error  may  be  estimated  from  the  spectrum  above  (uw,  v..)  for  the  detailed  grid  - 


-57- 


but  only  at  the  subset  of  the  "coarse"  frequencies.  For  practical  applications  it 
would  therefore  probably  be  more  valuable  simply  to  test  various  interpolation  pro¬ 
cedures  (e.g.  spline  functions)  for  Ah2-construction  and  the  subsequent  interpolation 
of  the  computed  far-zone  effects,  and  choose  the  method  with  least  high-frequency 
leakage.  Such  an  "optimum"  interpolation  method  is  necessary  anyway  in  order  to 
interpolate  results  from  the  FFT  computation  grids  to  actual  station  locations. 

Time  has  not  allowed  actual  implementation  of  the  FFT  methods  for  terrain 
effect  computations  within  the  present  project.  However,  recent  results  obtained 
by  Sideris  (1984)  seem  very  promising:  in  a  small  test  area  of  the  Rocky  Mountains 
(tc  range  4-22  mgal),  gravity  terrain  corrections  computed  with  FFT  showed  sub-mgal 
accuracies  when  compared  to  a  space-domain  prism  integration,  using  a  1km  x  1km 
elevation  grid. 

7.7  The  Linear  Approximation  and  Error  Studies 

In  addition  to  allowing  the  use  of  FFT  for  the  evaluation  of  general  terrain 
effects,  the  linear  approximation  also  comes  in  very  handy  in  the  study  of  error 
propagation,  e.g.  used  for  answering  questions  of  the  type:  given  a  certain 
statistical  behavior  of  the  gravity  field  and  the  topography,  to  what  extent  will 
it  be  beneficial  to  take  the  topography  into  account?  and  how  detailed  will  the 
height  information  be  needed?  etc. 

In  this  report  emphasis  is  on  residual  terrain  reductions  -  with  respect 
to  a  180  x  180  spherical  harmonic  expansion.  It  is  therefore  natural  to  work  with 
planar  (flat  earth  approximation)  error  analysis,  briefly  outlined  in  the  sequel 
as  it  is  not  too  familiar  to  many  geodesists.  The  basic  descriptor  of  the  statis¬ 
tical  properties  of  the  variations  of  the  gravity  field  (and  the  elevations)  is 
the  covariance  function,  e.g.  for  gravity  anomalies  at  a  reference  level. 


C(x,  y)  =  E{Ag(x',  y’)  Ag(x’+x,  y’  +  y)J 


(7.46) 
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where  E  is  the  ensemble  expectation  operator.  Within  the  area  of  interest  - 
approximated  by  an  infinite  plane  -  the  signal  Ug)  is  assumed  to  be  stationary. 
Fourier  transformation  of  C  yields  the  power  spectrum  (or,  rather,  power  spectral 
density) 

%,Ag(u’  v)  =  E  v)l2}  =  /  C(x’  y)el(ux+vy)dxdy  (7.47) 

For  an  isotropic  process,  C(x,  y)  =  C(r),  r  =  /F+y7,  (7.47)  will  as  earlier  men¬ 
tioned  be  a  Hankel  transform 

#  Ag  >Ag(  “)  =  2-tt  C(r )  -.2*  /“  rC(r)  JQ(tur)dr  (7.48) 

Of  special  importance  for  error  studies  is  Parsevals  formula  for  Hankel  transforms 
(Papoulis,  1968): 


r  r|  f ( r) |  2dr  =  J"  «|f((o)|2d«  (7.49) 

0  0 

Thus,  given  a  power  spectrum  <j»(io),  the  variance  of  the  signal  may  be  obtained  as: 

a2  =  C(Q)  =  /  a) <> ( cu)  dw)  (7.50) 

Spatial  extensions  are  obtained  by  upward  continuation  (6.14),  e.g.  for  the  potential 

d>TT(u),  zlt  z 2 )  =  <tTT(o) )e~w^zi+Z2'  (7.51) 

C-j-j(r,  Zj ,  z2)  =  -j— j-(o>  )e  Z2^Jo(ur)dr  (7.52) 


See  e.g.  (Nash  and  Jordan,  1978).  The  last  formula  -  which  is  simply  the  inverse 
transform  of  (7.48)  -  is  analogous  to  the  well-known  spherical  covariance  function 
expansion  in  Legendre  polynomials 

00  D  2  £+ 1 

CTTU,  r  ,  rj  =  j  o  [-r— ]  P,(cos*) 

1=2  i  rir2  1 


(7.53) 


where  is  the  degree  variances  ("spherical  power  spectrum"),  R  the  radius 
of  the  earth  reference  sphere  etc.  The  degree-variances  and  power  spectrum  are 
closely  related,  $md  a  unique  asymptotic  correspondence  exists  (Dorman  and  Lewis, 
1970) 


al  =  *TT( 


(7.54) 


More  details  will  be  given  in  a  subsequent  OSU  report. 

Power  spectra  and  covariance  functions  for  other  quantities -are  easily  derived 
by  using  the  expressions  of  the  quantities  in  frequency  domain,  e.g.  from  (6. 15 )- ( 6 . 17) 


%,Ag  S  “2*TT 


(7)  ♦ 


TT* 


nn 


PTT 


(7.55) 

(7.56) 


For  an  isotropic  field  (that  is  <t>Ag  A^  isotropic)  these  equations  and  Parsevals 


equation  gives  the  important  corollary  that 


aAg  =  'Y2(a!  +  =  2y2a|  (7.57) 

In  other  words,  the  gravity  variance  is  double  the  variance  of  each  of  the  de¬ 
flection  components.  Thus  0=  a ^  =  1"  corresponds  to  a Ag  -  6.7  mgal. 

To  describe  the  covariance  function  in  a  given  area,  simple  analytical  expres¬ 
sions  are  traditionally  used.  The  "Poisson"  and  "inverse  distance"  covariance 
functions  of  Moritz  ( 1980)  are  especially  important  for  gravity,  since  they  have  simple 
analytical  spatial  extensions.  For  topography,  however,  empirical  investigations 
of  U.S.  data  (cf.  next  section)  indicates  that  better  overall  fits  are  obtained 
using  exponential  covariance  functions  (so-called  first  order  Markov  models). 

The  basics  of  these  three  simple  models  may  be  outlined  as  (See  Moritz  (1980). 


Name 

'  C(r) 

0  (m) 

Correlation 

Length 

1)  "Poisson" 

a2 

[l+(r/D2)]3/2 

2irc2D2  e-a,° 

x,  =0.77  D 

2)  "Rec.  distance" 

a2  ! 

[l+(r/D)2]^ 

-«D 

2tto2D  - - 

CO  ! 

1 

x,  =1.73  D 

h 

3)  "Markov" 

e-r/D 

2^02D2  [1+a)2D2]3/2 

x,  =0.69  D 

These  functions  are  shown  in  Figure  15  together  with  an  actual  topographic  data 
example . 

7.8  Error  Studies  of  DTM  Resolution  Requirements 

To  give  an  example  of  error  analysis  in  terrain  reductions,  the  representation 
error  for  terrain  effects  on  gravity  and  deflections  at  altitude  will  first  be 
studied.  In  other  words,  the  resolution  requirements  for  a  digital  terrain  model 
to  give  terrain  effects  of  a  certain  accuracy  will  be  studied. 

Assume  topographic  mean  elevations  to  be  given  on  a  grid  of  spacing  (ax, Ay) . 

If  the  grid  elements  are  reasonably  "square",  then  the  mean  may  be  approximated 
with  a  mean  over  a  circle.  This  is  advantageous,  as  isotropy  then  will  be  "conserved". 
Then  to  first  order  the  terrain  effect  computed  from  the  mean  elevations  (neglecting 
known  station  elevations)  may  be  expressed  as: 

Ag'  =  2TrGP^hmean  *  2ttGp  J  Ah(x',  y')  dx'dy',  a  =  0.56/AxAy  (7.51) 

C 

where  Ah  =  h-h^  is  the  residual  elevation  and  C  a  circle  of  radius  a,  centered 

at  the  computation  point.  (7.51)  is  again  a  convolution 

5p  ( 1  r  <  a 

Ag'  =  —§■  (f*Ah),  f(r)  =  <  (7.52) 

0  (0  otherwise 


with  transfer  function  (Papoulis,  1968) 
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Figure  15  Empirical  power  spectrum  and  covariance ‘function  for  elevations  ah 
in  a  110x110  km  area  in  the  Smoky  Mountains,  eastern  USA,  and  fit  with  3  simpl 
two  parameter  covariance  models.  An  excellent  fit  is  obtained  with  the  expo¬ 
nential  (Markov  model ) . 


f(«)  =  2irf(u)  = 


(7.53) 


JL  is  the  Bessel!  function  of  order  1.  The  shape  of  the  function  f  is  shown 
in  Figure  16. 
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Figure  16 

By  upward  continuation  to  elevation  H  the  spectrum  then  becomes 


aV(u,  v)  =  4*Gp  ^4^-'  e‘wHAh(u,  v) 


(7.54) 


By  comparison  to  tne  "exact"  terrain  effect 


Ag(u,  v)  =  2irGp  e"uHAh(u,  v) 


(7.55) 


the  representation  error  e  =  Ag1  -  Ag  is  seen  to  be 


e(u,  v)  =  2irGp  (1  -  2  iliis*!)  e"wHAh(u,  v) 

do) 


(7.56) 


and  thus 


v<-> ■  (2’&>2  <■  - 


(7.57) 


For  deflections  of  the  vertical  essentially  the  same  formula  holds.  From  formula 
(7.21)  it  is  seen  that  the  deflection  error  expressions  corresponding  to  (7.56) 
will  be  simply 


^(u,  v) 

v> 


•  HD » ‘(u' #) 


(7.58) 
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For  an  isotropic  field,  <|>  =  .  we  get  directly  using  the  statistical  indepen¬ 

dence  of  c  and  n  ; 


<p£  e  M  =  0£  (“)  =  hl{\l  I2  +  |e  I2}  =  -^2  A  (u,  v) 

^  ^  n  n  ^  1  '  fce 


(7.59) 


Thus,  results  obtained  for  gravity  anomalies  in  the  error  analysis  may  be  directly 
applied  to  deflections  of  the  vertical  as  well,  when  the  previously  mentioned 
"conversion  factor"  of  6.7  mgal/arc  sec  is  used. 


To  get  error  estirates  an  exponential  ('Markov")  model  C  u  Ar)  =  02  e"r/D 

AhAh  Ah 

will  be  assumed  for  th<-  topography  covariance  function.  B/  (7.50)  and  (7.57) 
the  error  variance  wil  be 

a2  =  -J-  f  u  <p  U)  dw 

£  2lt  '  q  T££ 


(  2tt  Gp  ) 2  o2h  J“u(1 


Jllauil)2  -2.H  - Df - — -  du. 

[i+u2D2]3/2 


(7.60) 


It  is  convenient  to  normalize  the  parameters  with  respect  to  D,  introducing  a 

g  H 

dimensionless  averaging  parameter  a 1  -  -g-  and  elevation  H1  =  g.  Then  by  shift 
of  variable  t  =  wD: 


a2  =  (  2tt Gp  ) 2  a2h  F2  ( a  ’  ,  H1) 


(7.61) 


F2(a. ,  «■)  •  r  (i  -  2  =-2tH'  rAiw dt 

o  at  Ll+t  J ' 


(7.62) 


This  integral  has  no  simple  analytical  expression.  It  has  been  integrated  numerically, 
using  a  standard  adaptive  numerical  integration  subroutine  and  using  polynomial  ex¬ 
pansions  for  J;(t)  given  by  Abramowitz  and  Stegun  (1965).  The  result  (i.e.  the 
function  F,  square  root  of  7.62)  is  shown  in  Figure  17. 

The  r.m.s.  computation  error  op  at  elevation  H  should  ordinarily  be 
compared  to  the  "actual"  r.m.s.  terrain  effect  °  (H)  at  elevation  H.  The 


relative  r.m.s.  effect  at  altitude  will  be  given  simply  by  the  integral  (7.62) 
without  the  Bessel  function  term,  i.e.  as  the  limit  of  F  for  large  a'  , 
additionally  shown  on  Figure  17. 


Figure  17  R.m.s.  error  integral  (7.62)  for  terrain  reduction  of  gravity  and  deflec 
tions.  The  graphs  show  the  ratio  between  the  r.m.s.  computation  error  oe  and 
2irGp  aAh  as  a  function  of  normalized  averaging  radius  for  various  normalized  ele¬ 
vations.  Asymptotic  values  shown  at  left. 


To  give  an  example  of  application  of  the  error  curves,  consider  the  Smoky 
Mountains  area,  a  typical  "mild"  mountainous  area.  From  the  topography  covariance 
function  (Figure  15)  we  have  -  305  m,  correlation  length  -  7.6  km  and 
thus  D  =  11.0km.  First  consider  stations  at  the  level  of  the  topography,  H'  =  0. 
The  r.m.s.  variation  of  the  residual  terrain  effect  will  be  34  mgal  and  5.2"  for 
gravity  and  deflections  respectively.  To  compute  terrain  effects  with  a  6.6mgal/l" 
r.m.s.  error  (F  -  0.19),  Figure  16  gives  a'  -  0.008  and  by  (7.51)  ax  =  Ay  -  1.6  km 
Hence,  a  digital  terrain  model  with  a  grid  spacing  around  1.6  km  will  be  needed 
to  give  1"  -  deflections.  At  altitude,  say  H  =  10  km,  the  r.m.s.  terrain  effect 
is  seen  to  be  only  -13  mgal.  To  get  e.g.  a  1  mgal  error  (i.e.  F  -  1/34),  the 


figure  indicates  a'-  0.7  and  thus  the  resolution  of  the  digital  terrain  model 
needs  to  be  only  13-14  km. 

For  qeoid  undulations  similar  error  curves  may  be  computed.  In  the  linear 
approximation  height  anomaly  residual  terrain  effects  at  the  surface  of  the  topo¬ 
graphy  will  be 

c  =il&i4h  (7.63) 

Y  u> 

« 

Analogously  to  gravity  and  deflections  the  following  expressions  for  the  height 
anomaly  variance  and  error  variance  are  obtained 


o  _  /  2"n  Go  ^  2  o  p>9  c00  1 

a *  -  ( - )  of.  |  - 57- 

4  Y  Ah  J  0  a)[l+(1)2D2]3/2 


dw 


(7.64) 


2  =  (2jGp_\2  a2  D2  f°°  (1  -  2  J 1  ( aai ) ) 

e _  y  Ah  Jn  aw 


1 


)[l+w2D2]3/2 


dw 


(7.65) 


The  variance  a2  computed  by  (7.64)  is  infinite  (opposed  to  the  finite  c2  ). 

n 

This  is  a  phenomena  analogous  to  the  infinite  potential  effect  of  the  Bouguer 
plate:  very  long  wavelengths  in  the  topography  results  in  very  large  geoid  effects. 
However,  although  the  simple  Markov  covariance  model  used  has  energy  at  long  wave¬ 
lengths,  this  will  not  be  the  case  for  "real"  residual  terrain  reduction  with 

respect  to  a  180  x  180  reference  surface.  Ideally,  no  power  should  remain  below 

180 

the  "reference"  frequency  w0  =  — pp.  Thus,  a  better,  less  conservative  estimate 
of  the  "local"  height  anomaly  variation  is  obtained  by  integrating  from  w0  rather 
than  0  in  (7.64)  and  (7.65) 

a2  =  (2iG*L)2  a2  D2  g2  (D,  a')  (7.66) 

€ ^  Y  An  1 

G?(0,  (1  .  2  Jxiilti)2 

0 


where  the  substitution  t  =  <»)D,  a 


=  a/D,  has  been  used  again.  The  function 


Gt  is  shown  in  Figure  18.  Note  that  the  error  variance  is  nearly  insensitive 
to  the  omission  of  the  low  frequencies. 
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Figure  18  R.m.s.  error  integral  (7.67)  for  height  anomaly  residual  terrain  effects 
with  respect  to  a  180  x  180  spherical  harmonic  reference  surface. 

For  isostatic  reductions  the  "full"  height  anomaly  variance  will  be  well 
defined.  At  the  accuracy  levels  of  the  linear  approximation,  the  isostatic  reduc¬ 
tion  may  be  viewed  as  a  mass  plane  compensation  at  depth  T.  Then  the  isostatic 
reduction  transfer  function  will  be  simply 


(1  -  e"wT)  h 

y  u 


(7.68) 


where  it  is  now  the  elevation  h  and  not  the  residual  elevation  Ah  which  is 
used.  The  isostatic  error  variance  integral  then  becomes 


02  .  (*&.f  ^  „2  fi|  r) 


(a'.T')-r  (1  -  2  MlJAL)2  (1  .  e~tTV  —A 


(7.69) 


t[l+t2]3/2 


again  through  the  substitution  t  =  u>D,  T'  =  ^,  a'  =  The  dimensionless  function 
G2  is  shown  in  Figure  19. 


Figure  19  R.m.s.  error  integral  for  isostatic  height  anomaly  reductions.  Asymptotic 
values  for  large  a'  (i.e.,  the  variance  integral)  shown  in  brackets  at  right. 

For  a  RTM-reduction  example,  consider  again  the  Smoky  Mountains  area,  D  =  11  km, 

=  305  m.  From  Figure  18  and  (7.66)  the  r.m.s.  variation  of  the  RTM  height 

anomalies  is  seen  to  be  only  -  38  cm,  which  verifies  the  earlier  claimed  advantage 

of  the  RTM  reduction:  terrain  effects  on  the  geoid  are  very  small.  Now,  say  if  an 

accuracy  of  10  cm  is  wanted  for  the  geoid  terrain  effects.  Figure  18  shows  that 

a'  -  2  will  be  sufficient,  corresponding  to  a  necessary  gridspacing  of 
2 

&x  3  Ay  q — gg-  •  11  km  -  40  km,  again  demonstrating  the  insensitivity  of  geoid  terrain 
effects  to  the  very  local  topography. 

To  give  an  example  of  application  of  the  error  analysis  in  oceanic  areas, 
an  empirical  bathymetric  covariance  function  for  a  10°x  10°  trench  area  in  the 
Pacific  is  shown  in  Figure  20.  In  this  area  the  mean  depth  is  3.2  km,  o.  -  1553m 
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Figure  20  Power  spectrum  and  covariance  function  for  bathymetry  in  a  10°  x  10c 
area  around  the  Tonga  trench.  No  reference  elevations  have  been  subtracted. 
Again  the  exponential  Markov  model  is  seen  to  give  the  best  fit  of  the  simple 
covariance  models,  (power  spectrum  has  been  plotted  for  p  =  2.67.  To  get 
p  =  1.64  results,  simply  subtract  2  dB  from  given  values). 


and  Xjj  -  107  km,  giving  0  -155  km.  For  a  trench  area  the  isostatic  compensation 
is  known  to  be  partly  based  on  density  anomalies  in  the  upper  mantle,  and  therefore 
a  fairly  deep  isostatic  compensation  level  T  must  be  chosen.  Assume  e.g.  T=  40  km. 
Then  T'  -  0.26  and  from  Figure  19  and  (7.69)  with  p  =  -1.64  the  variance  of  the 
isostatic  geoid  effect  is'  seen  to  be  a  -  3.4  m  (a  direct  computation  of  the 


isostatic  effects  have  yielded  a  =  3.1  m,  thus  validating  the  error  analysis 
approach,  cf.  next  chapter).  To  compute  the  isostatic  geoid  effect  with  an  ac¬ 
curacy  comparable  to  SEASAT  satellite  altimetry,  say  a€  ~  10  cm,  then  we  must 
have  G2  -  0.006  and  hence  a1 -  0.09,  giving  necessary  grid  spacing  Ax-  25  km. 

An  "inverse"  example  may  also  be  given:  recently  a  global  5 ' x  5 ’  mean 
bathymetry  data  set  (SYNBAPS)  was  released,  covering  most  of  the  earths  oceans 
between  75°  N  and  75°  S.  In  the  "rough"  Tonga  trench  area,  the  SYNBAPS  data  cor¬ 
responds  to  a' -  0.033,  giving  an  r.m.s.  isostatic  computation  error  oe^~3  cm. 
However,  it  should  be  remembered  that  this  number  corresponds  to  the  linear  ap¬ 
proximation  and  mass  plane  condensation , and  assumes  the  mean  elevations  to  be 
error  free.  Therefore  SYNBAPS  derived  isostatic  geoids  might  be  significantly 
more  in  error,  naturally  especially  in  areas  of  poor  bathymetric  data  coverage. 


8.  Spectral  Characteristics  and  Covariance  Functions  for  Local  Topography  and 

Terrain  Effects 

In  the  present  section  key  parameters  describing  the  statistical  behavior 
of  the  local  topography  will  be  investigated  for  a  number  of  different  sample 
areas,  representing  various  types  of  topography,  from  nearly  flat  to  alpine  areas 
in  the  United  States. 

Power  spectrum  and  covariance  functions  have  been  estimated  from  available 
0 . 5 ' x  0.5'  elevations,  supplied  by  the  National  Geodetic  Survey,  using  a  simple 
FFT  approach. 

With  this  approach,  the  residual  topography  power  spectrum  (power  spectral 
density)  is  obtained  from  the  discrete  fourier  transform  (7.37)  of  Ah  by 


where  Ax,  Ay  are  grid  spacing  of  the  nxm  given  elevation  grid.  The  power 
spectrum  ^  will  be  defined  in  a  frequency  square  between  the  Nyquist  fre- 


quencies  (7.40)  tujj,  ±Vfl.  It  will  be  symmetric  with  respect  to  the  zero-frequency: 
4>(u,  v)  =  4>(-u,  -v).  By  the  inverse  fourier  transform  (7.38)  the  2-dimensional 
covariance  function  is  obtained: 


CAhAh^x>  y)  =  ^nm^  ^AhAh^u* 


(8.2) 


This  function  will  be  defined  within  the  square  ±^  ax,  Ay,  and  will  again  exhibit 
symmetry  with  respect  to  (0,0). 


In  the  sequel  <j>  and  C  will  be  investigated  in  a  number  of  cases  for  both 
topography,  terrain  corrections  and  observed  gravity  and  geoid.  Since  we  are 
mainly  concerned  with  isotropic  processes,  the  results  will  be  averaged  along 
circles,  to  give  the  "isotropic"  covariance  function  C(r)  and  power  spectrum 
<j>(oo).  However,  to  get  an  idea  of  the  anisotropy,  a  contour  plot  of  the  2-dimen¬ 
sional  covariance  function  will  be  given  as  well.  For  gravity  also  degree-vari¬ 
ances  will  be  given,  based  on  formula  (7.54). 

The  power  estimate  (8.1)  may  be  improved  by  taking  the  finite  extent  of  the 
given  data  into  account  through  the  use  of  window  filters.  For  tests,  a  two- 
dimensional  Hanning  window  has  been  applied  in  the  frequency  domain  (see  e.g. 
Engelis,  1983),  the  main  effect  being  a  less  "noisy"  <p(u,  v).  The  discrete  Hanning 
window  in  the  frequency  domain  is  nothing  but  a  simple  smoothing,  with  a  filtered 
frequency  value  given  as  a  weighted  mean  of  the  value  itself  and  the  adjacent 
frequencies.  Since  this  smoothing  already  to  some  degree  is  performed  by  the 
radial  smoothing  (and  the  prime  interest  is  in  the  overall  shape  of  the  curves) 

I  have  found  no  practical  need  for  windowing,  and  thus  the  spectral  estimates 
represent  "raw",  unfiltered  FFT-results. 


Examples  of  results  have  already  been  shown  in  Figures  15  and  20.  The  power 
spectra  will  be  given  in  units  of  mgalzdegree2  (1  degree-  111  km),  where  elevations 
have  been  converted  to  gravity  using  the  simple  Bouguer  factor  for  p  =  2.67  g/cm3. 


(8.2) 


gm  =  2irGpAh  -  0.1119  [mgal/m]  ■  Ah 
and  results  are  shown  logarithmic  in  dB 

<j>  [dB]  =  10  log1Q  U  [mgal 2degree2])  (8.3) 

as  a  function  of  the  frequency  v  =  -!j~  in  units  cycles/degree. 

The  anisotropy  will  be  indicated  through  a  small  contour  plot  of  the  normalized 
2-dimensional  covariance  function  at  levels  C(r)  =  1.0,  0.8,  0.6,  ...  (Figure  21). 

To  get  a  quantative  measure  of  the  anisotropy,  an  anisotropy  index  will  be  defined  as 

maximal  Xj, 

anisotropy  index  -  min1mal  (8.4) 

where  X,  is  the  correlation  length  in  a  certain  direction,  i.e.  the  distance 
from  origo  for  which  the  covariance  is  half  its  zero-value  CQ.  An  isotropic 
field  will  have  an  anisotropy  index  of  1. 

The  second-order  gradient  variance  GQ  will  be  used  in  the  sequel  as  a  measure 
of  high-frequency  content.  As  shown  by  Moritz  (1980),  Gq  will  be  given  as  the 
curvature  of  the  (gravity)  covariance  function  at  origo.  It  will  be  estimated 
directly  from  the  empirical  covariance  function,  using  a  symmetrical  spline  inter¬ 
polation  procedure.  The  result  must  be  used  with  some  caution,  especially  for 
the  topography  where  the  applied  Bouguer-approximation  is  very  crude  for  second-order 
gradients.  Since  the  spline  curvature  determination  will  be  directly  dependent 
on  the  data  grid  spacing,  G Q  will  merely  represent  a  lower  limit  of  the  gradient 
variance,  giving  the  gradient  variance  of  some  smooth,  filtered  version  of  the 
actual  gradient  field.  The  gradient  variances  will  be  given  in  E2, 

IE  =  10"9s2  =0.1  mgal/km. 
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Figure  21  Empirical  residual  elevation  covariance  function  for  the  Colorado 
Rocky  Mountain  test  area  of  Figure  12  (l°xl°  block,  latitude  39°  to  40°  N, 
longitude  107°  to  106°  W). 

TOP:  2-dimensional  covariance  function  from  FFT, 

RIGHT:  same,  but  shown  with  0.2  C0  contour  interval  (and  only  to  ±15  arcmin), 
BOTTOM:  radially  averaged  covariance  function.  Poisson  and  reciprocal  distance 
covariance  models  also  shown  (again  an  exponential  model  would  give  a  better 
fit,  cf.  Figure  15) 
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8.1  Topographic  Covariance  Functions  for  U.S.  Sample  Areas 

A  number  of  l°xl.2°  nearly  square  blocks  (each  containing  120x144  elevation 
points)  have  been  analyzed  using  0.5'x0.5'  (-1  km)  elevations  from  the  NGS  data 
base.  Each  analysis,  with  FFT,  plotting  etc.,  used  c.5  sec  CPU-time  on  The  Ohio 
State  University's  Amdahl  computer. 

The  areas  are  shown  in  Figure  22,  and  represent  various  types  of  topography: 
alpine  (Colorado,  Sierra  Nevada),  alpine- mountainous  with  volcanoes  (Washington), 
mountainous  plains  (New  Mexico),  older  mountains  (Smoky  Mountains)  and  flat  to 
hilly  lowland  (Ohio).  The  two  1°  New  Mexico  blocks  comprise  most  of  the  "White 
Sands"  New  Mexico  test  area  for  comparison  of  gravity  field  modelling  techniques, 
see  e.g.  (Tscherning  and  Forsberg,  1983). 

All  elevations  analyzed  are  residual  elevations  Ah  =  h  -  h  where  the 
180x180  spherical  harmonic  expansion  of  the  earth's  topography  of  Rapp  (1982) 
has  been  used  as  reference  elevation.  The  Bouguer-deri ved  topographic  gravity 
is  thus  to  first  order  the  RTM  180  terrain  effect. 

Results  are  shown  in  Figures  23,  24,  and  Table  4.  Despite  the  different 
types  of  topography,  the  results  are  amazingly  similar,  with  the  prime  variation 
being  in  the  r.m.s.  variation  of  CQ,  less  in  the  shape  of  the  curves  and  the 
correlation  length.  The  correlation  length  varies  in  the  range  6.5-12.2  km,  while 
the  r.m.s.  residual  elevations  varies  from  45  m  to  785  m.  For  the  eight  test 
areas,  maximal  and  minimal  roughness  of  RTM-terrain  effect  field  is  obtained  in 
the  Sierra  Nevada  and  Ohio  blocks  respectively:  for  gravity  the  r.m.s.  terrain 
effect  will  be  88  and  5  mgal  respectively,  for  deflections  13"  and  .8"  and  for 
height  anomalies  96  cm  and  5  cm,  using  the  results  of  section  7.8.  The  degree 
of  anisotropy  is  seen  in  some  cases  to  be  nil,  in  other  cases  rather  severe. 

Since  the  observed  gravity  field  will  be  dominated  by  the  topographic  effects 
in  mountainous  blocks,  the  same  degree  of  anisotropy  will  be  expected  in  the 
gravity  covariance  function,  thus  stressing  the  need  for  utilization  of  terrain 
reduction  methods. 
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Figure  22  l°xl.2°  sample  areas  for  topographic  covariance  functions  and  power 
spectra  (hatched  small  rectangles)  and  gravity  analysis  sample  areas  (large 
rectangles)  in  the  United  States. 
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Figure  23  Empirical  power  spectra  (top)  and  covariance  functions  (bottom  and 
right)  for  residual  topography  in  selected  l°xl.2°  mountainous  blocks  of  western 


Table  4 

Statistical  Parameters  for  Residual  Topography,  l°xl.2°  Blocks,  0.5'x0.5'  Sampling 


39-40,  £3.1-81.9 
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8.2  Magnitude  of  the  Gravimetric  Terrain  Correction  -  Colorado  Area 


In  the  previous  sections  emphasis  has  been  on  the  dominant  terrain  effect 
on  gravity  anomalies  -  the  "Bouguer"  term  -2irGpAh.  To  get  the  complete  terrain 
effect,  the  additional  terrain  correction  is  needed.  However,  it  will  usually 
be  one  order  of  magnitude  smaller  than  the  Bouguer  term  -  in  the  present  section 
the  actual  magnitudes  will  be  investigated  for  the  Colorado  areas. 
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Figure  25  Terrain  corrections  (top)  and  free-air  anomalies(  bottom)  for  gravity 
stations  in  a  Colorado  area  (lat.  38.5-39°  N,  Ion.  106.5-106°  W).  Note 

"e  different  scale  used  in  the  plots. 
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Figure  25  shows  an  example  of  actual  gravity  data  and  terrain  corrections 
for  a  small  (h°xh°)  block  (NGS  data).  The  terrain  corrections  are  seen  to  be 
in  the  range  of  0-40  mgal,  opposed  to  a  Bouguer  RTM-effect  range  of  roughly  -120 
to  120  mgal.  For  such  a  rather  small  area,  the  free-air  anomalies  correlate  well 
with  elevation,  the  slope  being  essentially  a  measure  of  the  topographic  density. 
However,  since  the  terrain  corrections  have  a  tendency  to  be  large  for  the  higher 
stations,  the  free-air  slope  will  give  too  low  density  values:  in  the  case  of 
Figure  25  p-  2.1,  which  is  unrealistic.  Therefore,  while  terrain  corrections 
might  be  neglected  in  connection  with  e.g.  error  studies,  omission  in  other  cases 
implies  serious  systematic  errors,  especially  when  used  for  density  determination 
as  outlined  in  Section  5. 

To  study  further  the  terrain  correction,  computation  of  terrain  corrections 
was  done  in  a  2 ' x  3 ‘  grid  for  the  two  l°x  1°  Colorado  blocks  of  the  last  section 
(and  on  a  2 ' x  2 '  grid  for  the  h°xH°  area  of  Figure  25),  using  the  prism  integration 
programme  of  the  appendix.  The  choice  of  the  rather  coarse  computation  grid  size 
was  necessary  to  avoid  excessive  computation  times  (on  the  OSU  Amdahl  system 
CPU-time  requirements  were  c.  0.2  sec/station).  Figure  26  shows  the  result  for 
the  northeastern  block.  By  comparing  to  the  topographic  map  (Figure  12)  it  is 
clear  that  some  degree  of  undersampling  is  done  by  only  computing  on  a  2 ‘ x  3 '  grid, 
so  some  aliasing  might  be  expected  in  the  power  spectra  of  the  terrain  correction 
signals,  shown  in  Figure  27. 

At  the  2'  resolution  level,  the  terrain  correction  covariance  functions  (Figure 
27)  might  be  described  as  a  "white  noise"  signal  with  very  short  correlation  length 
plus  a  smaller  long  wavelength  signal  (remote  zone  effects).  Compared  to  the 
main  "Bouguer"  terrain  effect,  the  terrain  corrections  are  indeed  seen  to  be  smaller, 
as  shown  in  Table  5  -  however,  large  outliers  and  non-normal  error  distribution 
makes  it  important  to  use  the  best  possible  terrain  corrections  for  gravity  field 
modelling,  using  even  more  detailed  elevation  data  than  the  presently  used  0.5'x0.5' 
heights . 


Figure  26  Terrain  corrections  in  mgal  for  "Colorado  I"  l°x  1°  block,  computed 
on  a  2 1 x  3 '  grid.  A  topographic  map  of  the  same  area  is  shown  in  Figure  12. 
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Figure  27  Power  spectra  and  covariance  function  for  the  gravimetric  terrain 
correction  in  two  l°x  1°  blocks  of  Colorado.  The  main  terrain  effect  additionally 
shown  with  dots,  cf.  Figure  23. 
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Area 

RTM  Bouguer  term 
(mgal ) 

Y 

r.m.s.  h 

Terrain  corrections 

mean  std.  dev.  max 

[mgal ) 

I 

(39-40°  N,  107-106°  N) 

47.7 

3.5' 

6.4 

4.4 

55 

1.8' 

II 

(38-39°  N,  108-107°  N) 

48.5 

6.1' 

4.8 

3.7 

29 

2.4' 

III 

(38.5-39°  N,  106.5-106°  W) 

51.9 

4. 6  ‘ 

6.4 

4.4 

27 

2.4' 

Table  5  Gravity  terrain  corrections,  Colorado  areas. 

The  obtained  mean  terrain  corrections  may  be  used  to  test  Siinkel's  approx¬ 
imate  formula  (cf.  Section  7.2): 

Ou 

mean  terrain  correction  >  3itGp  (8.5) 

For  the  three  areas,  values  of  4.7  mgal ,  2.8  mgal  and  4.2  mgal  are  obtained. 

The  mean  terrain  corrections  obtained  by  (8.5)  are  thus  seen  to  be  somewhat  too 
low  (this  deficiency  is  primarily  a  consequence  of  the  tendency  of  the  topographic 
covariance  functions  to  be  exponential:  for  such  covariance  function  the  integral 
formula  underlying  (8.5)  is  undefined,  cf.  Siinkel ,  1981a,  p.  62). 

8.3  RTM  Geoid  Effects  -  Colorado  Area 

When  using  a  residual  terrain  reduction  with  respect  to  a  180x180  spherical 
harmonic  elevation  reference  surface,  the  study  of  section  7.8  showed  the  terrain 
geoid  effects  to  be  fairly  small.  In  this  section,  an  example  of  actual  magnitudes 
encountered  in  an  extreme  case  -  the  mountainous  part  of  Colorado  -  will  be  given. 
The  4°x  4°  block  of  Figure  22  -  comprising  essentially  all  of  Colorado  west  of 
Denver  -  will  be  studied  based  on  4 ' x  5 *  mean  elevations. 

Figure  28  shows  the  residual  topography  covariance  function  for  the  4'x5' 
elevations.  RTM  geoid  undulations  were  computed  from  this  grid  in  a  "fixed-area" 
reduction  taking  only  into  account  the  4°x  4°  square.  The  geoid  effects  varied 
from  1  to  -3  m  with  a  mean  of  -80  cm  and  an  r.m.s.  variation  of  87  cm  with  respect 
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Figure  28  Covariance  functions  for  residual  topography  (lower)  and  computed  RTM-180 
geoid  effects  (upper)  for  the  mountainous  part  of  Colorado.  The  180x180  reference 
elevation  surface  is  seen  to  give  a  good  fit  in  Colorado,  the  height  covariances 
being  near  zero  for  ^>1°. 

to  the  mean  (Figure  29),  corresponding  to  indirect  effects  on  gravity  anomalies 
below  1  mgal.  For  comparison,  using  the  statistical  study  of  Section  7.8,  the  ele¬ 
vation  covariance  parameters  (aAh  =  369  m,  =  8.4')  yields  an  r.m.s.  geoid 
variation  of  66  cm,  the  slightly  too  low  value  being  primarily  due  to  the  remaining 


small  long-wavelength  geoid  effects,  evident  from  Figure  28. 


9.  Terrain  Reductions:  Spectral  Characteristics  and  Covariance  Functions  for 
Local  Gravity 

By  using  known  density  anomalies  in  gravity  field  modelling,  the  purpose 
is  to  obtain  a  more  smooth  field,  more  suitable  for  interpolation  and  prediction. 

In  this  section  the  actual  smoothing  obtained  using  the  residual  terrain  reduction 
will  be  studied  for  a  number  of  U.S.  sample  areas,  shown  in  Figure  22,  each  repre¬ 
senting  various  types  of  topography  and  geologic  setting:  the  2°x  2°  Sierra  Nevada 
block  contains  the  highest  part  of  the  Sierra  Nevada  mountains  plus  a  part  of 
the  California  Valley,  and  is  characterized  by  a  large  horizontal  gravity  gradient, 
relating  to  changes  in  crustal  and  upper  mantle  structures.  The  4°x  4°  Colorado 
block  has  high  mountains  all  through,  and  a  thick  crust  giving  rise  to  very  low  Bouguer 
anomalies.  The  4°x  4°  New  Mexico  block  has  both  mountains  and  plateau-type  topo¬ 
graphy,  and  major  density  anomalies  associated  with  the  Rio  Grande  rift  system. 
Finally,  the  4°x  4°  Ohio  area  represents  lowlands  and  an  area  geology  of  primarily 
thick  paleozoic  sediments. 

The  4  areas  have  a  reasonable  coverage  of  terrestrial  gravity  anomalies  in 
the  NGS  data  base.  To  obtain  covariance  functions  and  power  spectra  of  the  data, 
the  terrain-corrected  Bouguer  anomalies  were  gridded  in  a  4'x  5'  (ca.  7. 5x7. 5  km) 
grid,  using  a  truncated  collocation  gridding  algorithm  (Cruz,  1983)  where  the 
value  at  a  point  is  obtained  from  the  5  closest  neighboring  points  (used  subsequent 
to  an  initial  thin-out  screening,  where  only  one  data  point  per  l'xr  "pixel" 
was  retained).  An  example  of  the  obtained  Bouguer  anomaly  grid  is  shown  in  Figure  30 
To  obtain  RTM-gravity  anomalies,  relating  to  a  180x180  spherical  harmonic 
expansion,  use  is  made  of  (7.5): 

RTM  anomaly  Agc  =  &g  -  AgRTM180  *  Ag  -  2ttGo  (h  -  hrgf)  +  tc  (9.1) 

For  local  gravity  field  modelling,  a  180x180  spherical  harmonic  gravity  reference 
field  corresponds  as  earlier  mentioned  to  the  RTM  reduction.  Thus,  by  (9.1)  the 
terrain-reduced  gravity  anomalies,  referred  to  a  180x180  reference  gravity  field. 


is  simply  obtained  by 


,reSldual  -  <49  -  49„f)  -  49C  -  49BOU9',er  -  49B^9Uer 


(9.2) 


i.e.  as  the  difference  between  the  usual  terrain-corrected  Bouguer  anomalies  and 
the  reference  field  Bouguer  anomaly  Ag^guer=  Agrgf  -  2-rrGphrgf,  cf.  Figure  31. 

From  the  gridded,  residual  gravity  anomalies,  gridded  free-air  anomalies  have 
been  reconstructed  by  a  simple  "inverse"  Bouguer  reduction,  neglecting  the  terrain 
correction,  which  -  as  demonstrated  in  section  8.2  -  is  small  compared  to  the 
"main"  terrain  effect. 

The  data  have  been  spectral  analyzed  using  FFT,  to  obtain  power  spectra  and 
covariance  functions,  as  described  in  section  8.  The  power  spectra  have  additionally 
been  converted  to  normalized  degree-variances,  using  the  formula  (7.54).  This 
allows  a  comparison  of  the  variability  of  the  gravity  field  to  global  spherical 
harmonic  degree-variance  models,  e.g. 


Kaula's  rule:  a.  -  0.7  •  10“10 


(9.3) 


Tscherning-Rapp' s  model:  -  4.47  •  10"10 


(0.999617)*+1  (9.4) 


see  e.g.  Moritz  (1980).  These  models  may  be  viewed  as  "average"  earth  models, 
and  for  a  particular  area  they  may  be  used  as  indications  of  the  smoothness  of  the 
gravity  field,  compared  to  the  global  behavior. 

An  example  of  the  obtained  two-dimensional  covariance  function  is  shown  in 
Figure  32,  results  are  shown  for  Colorado  in  Figure  33  and  for  the  other  areas 
in  Figure  34.  Table  6  outlines  the  main  statistical  parameters,  analogous  to 
Table  4. 

The  results  are  seen  to  confirm  what  intuitively  should  be  expected:  in 
mountainous  areas  the  use  of  terrain  reductions  significantly  reduces  the  variance 
of  the  gravity  field,  especially  for  the  shorter  wavelengths,  where  the  decrease 


Figure  31  Spherical  harmonic  reference  Bouguer  anomalies  in  Colorado,  obtained 
from  gravity  and  topography  expansions  of  Rapp  (1981,  1982),  complete  to  degree 
and  order  180.  Contour  interval  5  mgal. 


Figure  32  Two-dimensional  covariance  function  for  the  Colorado  block  for  terrain- 
reduced,  residual  gravity  anomalies  (argument  interval:  -2°  to  2°  in  latitude 
and  longitude). 

in  power  in  the  present  test  areas  is  around  12  dB,  corresponding  to  a  factor 
of  c.16  in  power  or  4  in  r.m.s.  variation.  After  reduction  the  spectra  are  fairly 
similar  to  the  lowland  gravity  spectrum  of  Ohio,  although  they  still  contain  some¬ 
what  more  energy  (this  should  not  be  too  surprising,  since  mountainous  areas  naturally 
are  areas  of  tectonic  activity  and  thus  large  geologic  density  anomalies). 

Also,  the  correlation  length  is  seen  to  increase  and  the  anisotropy  seems 
to  decrease  significantly  for  terrain-reduced  data,  quite  as  expected.  A  typical 
correlation  length  for  the  RTM-reduced  data  seems  to  be  around  15'. 

Comparing  to  the  globcl  degree-variance  models  (which  in  principle  also  con¬ 
tains  the  effect  of  an  "average"  topography)  the  reduced  gravity  data  lies  sig¬ 
nificantly  below  in  power  for  the  higher  wavelengths.  This  illustrates  the  general 
need  for  always  "adapting"  covariance  functions  to  particular  gravity  field  modelling 
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Figure  34  Degree  variances,  power  spectra  and  covariance  functions  for  gravity 
data  for  U.S,  test  areas.  (Note:  scale  of  anisotropy  plots  A  and  B  are  doub'p 
the  scale  of  C,  0,  and  E.) 
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area:  the  global  covariance  models  must  generally  be  modified  to  have  less  power 
at  higher  frequencies.  For  the  free-air  gravity  data,  only  the  New  Mexico  block 
is  seen  to  be  of  fairly  typical  "global"  type.  This  block  is  thus  a  good  repre¬ 
sentative  test  area  for  gravity  field  modelling  experiments. 

The  general  conclusions  of  this  section  supports  the  results  of  the  extensive 
study  of  North  American  gravity  covariance  functions  by  lachapelle,  Mainville 
and  Schwarz  (1983).  That  study,  based  on  5 ' x  5 '  gravity  data  and  using  a  completely 
different  approach  than  here,  shows  parameters  in  good  agreement  with  relevant 
parts  of  Table  6.,  e.g.  for  "block  11"  of  that  study,  covering  most  of  Colorado,  a 
free-air  CQ-value  of  2071mgal2is  given,  opposed  to  the  2061  mgal2  obtained  here. 

For  the  gravity  gradient  variance  GQ,  the  average  for  "block  11"  is  c.  1750  E2 
compared  to  4620  E2  of  Table  6  for  unreduced  data,  and  c.  400  E2  versus  324  E2 
for  "terrain-reduced"  data. 

These  GQ- values  are  probably  much  too  low,  the  spacing  of  the  data  being 
simply  too  large.  An  example:  for  the  2°x  2°  Colorado  area  38-40°  N,  107-105°  N 
(which  has  a  dense  gravity  coverage),  a  computation  of  Gq  based  on  2'x2.5'  data 
have  yielded  GQ-values  of  1410°  and  1006  E2  for  unreduced  and  reduced  data  respec¬ 
tively  -  and  by  utilizing  the  topography  results  of  Table  4,  a  0.5'x0.5‘  data 
grid  would  be  predicted  to  yield  GQ-values  around  60000  E2  for  free-air  anomalies! 
This  shows,  that  the  GQ-quantity  is  nearly  meaningless  in  mountainous  areas  - 
only  combined  with  a  suitable  filtering  through  well-defined  mean  value  operators 
or  upward  continuation  do  the  second-order  gradients  have  well-behaved  statistical 
parameters. 


10.  Isostatic  Reductions  of  Satellite  Altimetry  Data  In  Two  Areas  of  the  Pacific 

Many  of  the  features  of  the  ocean  bathymetry  are  of  such  dimensions,  that 
use  of  a  kind  of  "residual"  terrain  effect  with  respect  to  a  180x180  reference 
surface  would  simply  miss  the  bulk  of  the  geoid  effects  associated  with  these 
features.  In  this  section,  emphasis  will  therefore  be  on  isostatic  reductions 
and  their  relationship  to  geoid  heights  determined  by  satellite  altimetry. 

A  recently  released  global  bathymetric  data  set  -  SYNBAPS  -  covering  nearly 
all  of  the  earth's  oceans  from  75°  N  to  75°  S  with  5 * x  5 1  mean  depths,  provides 
a  very  convenient  data  set  for  the  computations  of  ocean  isostatic  geoid  effects. 

As  investigated  in  section  7.8,  the  resolution  of  the  SYNBAPS  data  allows  ocean 
"terrain  reductions"  to  be  computed  with  r.m.s.  errors  at  the  cm-level.  An  example 
of  the  SYNBAPS  data  is  shown  in  Figure  35. 
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Figure  35  SYNBAPS  5 ' x  5 ‘  Bathymetry  in  the  10°x  10°  Tonga  Block 
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The  bathymetric  data  was  used  to  compute  isostatic  geoid  undulations  by  the 
"tc"-programme  system  for  two  10°x  10°  areas  of  the  Pacific: 

"TONGA":  a  trench  area  of  extreme  bathymetry,  shown  in  Figures  35  and  36. 
and 

"TAHITI":  a  mid-plate  island  area  to  the  NE  of  Tahiti,  shown  in  Figure  37. 

For  both  these  types  of  areas  traditional  Airy  isostasy  is  expected  to  have 
limited  applicability  -  especially  for  the  trenches,  where  isostatic  compensation 
is  partially  absent,  and  the  existing  compensating  density  anomalies  occur  at 
rather  deep  levels,  associated  with  downgoing  slabs  etc.  The  failure  of  "traditional 
isostasy  at  trenches  is  evident  from  a  simple  numerical  consideration:  for 
a  normal  crust  thickness  T  =  30  km  and  density  contrast  of  0.6  g/cm3,  commonly 
applied  continental  values,  the  isostatic  mantle  "anti-roots"  will  end  above  the 
ocean  bottom  at  depths  greater  than  8000  m!  Also,  the  implied  average  oceanic 
crustal  thickness  of  13-14  km  is  too  thick  by  a  factor  of  nearly  2. 

Isostatic  reductions  are,  however,  very  useful  as  an  empirical  mathematical 
tool,  irrespective  of  whether  or  not  the  physical  reality  may  be  described  by 
simple  models.  In  the  sequel,  three  simple  types  of  isostatic  compensation  will 
be  tested: 

A)  conventional  Airy  isostasy,  T  =  30  km,  Ap  =  0.6  (With  average  oceanic 
depths  around  4000  m,  this  corresponds  to  some  degree  to  a  mass  plane 
condensation  T -  12  km). 

B)  a  compensation  in  the  upper  mantle  in  depth  ranges  20-60  km,  described 
as  a  mass  plane  compensation  at  depth  T  *  40  km. 

C)  compensation  in  deeper  parts  of  the  upper  mantle  and  lithosphere,  ap¬ 
proximated  by  a  mass  plane  compensation  at  depth  T  =  70  km. 

For  each  of  the  areas,  bathymetry  will  be  taken  into  account  for  a  fixed  area 
extending  2°  outside  the  10°x  10°  area,  i.e.  for  a  14°x  14°  area.  By  using  such 
a  fixed-area  isostatic  reduction,  the  computed  effects  will  have  a  varying  (arbi¬ 
trary)  bias.  It  is  therefore  only  the  variation  with  respect  to  the  mean  which 
is  indicative  of  the  "fit"  of  the  isostatic  models. 
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Figure  36  15 ' x  15 ’  bathymetry  in  the  Tonga  block.  Contour  interval:  500  m. 

Geoid  unaulations  for  the  two  areas  were  derived  from  the  combined  SEASAT-GEOS-3 
satellite  altimeter  data  base,  available  at  The  Ohio  State  University  (Liang, 

1983).  Sea  surface  heights  with  respect  to  GRS80  were  identified  with  geoid 
undulations,  and  weregridded  in  a  15 ‘ x  15 '  grid  using  the  truncated  collocation 
procedure  (5  closest  points),  analagous  to  the  gravity  gridding  of  section  9. 


Figure  37  1 5 ' x  15 '  bathymetry  in  the  Tahiti  block  (Tahiti  is  in  lower  left  and 
the  Tuamotu  Islands  are  in  the  center).  Contour  interval:  500  m. 

The  obtained  altimeter  geoids  are  shown  in  Figures  38  and  39.  The  correlation 
to  the  bathymetry  is  obvious  in  both  cases,  the  Tonga  Trench  having  e.g.  an  impres¬ 
sive  asymmetric  negative  geoid  anomaly  of  nearly  20  m,  whereas  the  islands  have 
associated  positive  geoid  anomalies  of  some  meters.  Degree  variances,  power  spectra 


Figure  38  15 ' x  15 '  combined  SEASAT/GEOS-3  satellite  altimetry  geoid  in  the  Tonga 
area,  referred  to  GRS80.  Contour  interval  0.5  m.  Altimeter  sea  surface  height 
points  used  in  prediction  of  gridded  values  shown  with  dots. 

and  covariance  function  for  the  geoids  are  shown  later  in  Figure  47. 

Topographic/isostatic  effects  computed  from  the  "shallow"  Airy  isostasy  (A), 
the  intermediate  depth  compensation  (B)  and  deep  compensation  (C)  for  the  Tonga 
area  are  shown  in  Figures  40-42.  From  these  figures  it  is  seen  that  the  trench 
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Figure  39  SEASAT/GEOS-3  altimeter  geoid  in  the  Tahiti  block 

itself  is  partially  uncompensated,  i.e.  a  very  deep  compensation  level  is  needed 
to  reproduce  the  large  geoid  anomaly  across  the  trench.  This  is  also  shown  in 
Figure  43.  On  the  other  hand,  some  features  (e.g.  most  of  the  islands)  seem  to 
be  compensated  at  much  more  shallow  levels.  There  is  therefore  no  simple  general 
isostatic  principle  which  is  applicable  in  this  area. 
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Figure  40  Isostatic  g eoid  "A"  in  the  Tonga  area,  conventional  Airy  isostasy  (power 
spectra  and  degree  variances  shown  in  Figure  47) 

For  the  Tahiti  block,  an  example  of  an  isostatic  geoid,  derived  from  the  bathymetry 
is  shown  in  Figure  44.  By  comparing  to  the  satellite  gecid  (Figure  39),  it  is 
seen  that  the  central  Tuamoto  Islands  area  is  well  represented  by  the  conventional 


Figure  41  Isostatic  geoid  "B"  in  the  Tonga  area  (mass  plane  compensation  in  depth 
40  km) 

Airy  isostasy,  while  the  smaller  Tahiti  islands  apparently  are  compensated  at 
deeper  levels  (Figure  45),  or  -  rather  -  is  partially  uncompensated.  By  subtracting 
the  computed  isostatic  geoids  from  the  altimetry  geoids,  smoother  "terrain-reu^ced" 
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Figure  42  Isostatic  geoid  C  in  the  Tonga  area  (mass  plane  compensation  in  depth 
70  km) 

"  Figure  43  Interpolated  observed 

sea  surface  heights  and  computed 
isostatic  geoid  profiles  across  the 
Tonga  trench  in  a  W-E  profile  at 
18°  S  (arbitrary  bias). 
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Figure  44  Airy  isostatic  geoid  ("A")  in  the  Tahiti  block.  (Power  spectrum  and 
degree  variances  shown  in  figure  47). 


\  Figure  45  Satellite  altimetry  geoid 
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Table  7  Geoid  undulation  variation  for  various  isostatic  reductions. 


(A:  T  - 

12  km,  B:  T 

=  40  km 

,  C:  T  = 

70  km,  massplane  comp 

) 

Geoid 

Variation 

Isostatic  Geoid 
std.  dev. 

Reduced  Geoid 
std.  dev. 

Area 

Geoid 

mean 

std.  dev. 

ABC 

A 

B 

C 

Altimetry 

37.61 

9.34 

7.93 

6.59 

4.84 

TONGA 

AU.-GEM9 

-0.18 

3.60 

1.61  3.12  5.17 

2.91 

2.97 

4.21 

Altimetry 

-2.80 

3.15 

3.41 

4.00 

4.82 

TAHITI 

AU.-GEM9 

-0.19 

1.69 

1.18  2.16  3.49 

1.04 

1.23 

2.28 

By  inspecting  the  numbers  of  the  table  it  is  seen  that  the  "smoothing"  effect 
of  the  terrain  reduction  is  not  very  marked,  and  it  is  apparently  mainly  related 
to  the  shorter  wavelengths  (alt.-GEM9).  However,  by  inspecting  the  isostatic 
geoids  and  comparing  them  to  the  "actual"  geoid,  it  is  clear  that  many  features 
of  the  ocean  geoid  are  essentially  nothing  but  bathymetric/isostatic  effects,  but  they 
do  not  show  up  very  well  in  the  statistics  because  the  geoid  signal  is  dominated 
by  effects  from  deeper  sources  (cf.  Figure  47).  Also,  error  sources  such  as  orbit 
errors,  sea  surface  topography  and  errors  in  SYNBAPS  might  play  some  role. 

From  the  statistics  there  seems  to  be  a  slight  favorization  of  the  conventional 
Airy  isostasy,  in  spite  of  the  obvious  deficiencies  for  the  trench  (this  would 
be  a  good  case  for  utilizing  a  "geologic"  density  anomaly  model,  e.g.  a  dipping 
slab).  For  global  studies,  the  geoid  anomalies  observed  over  trenches  carry  a 
very  high  weight  due  to  their  magnitude,  and  when  solving  for  "optimal"  compen¬ 
sation  depths  somewhat  too  deep  levels  might  be  obtained..  A  recent 
global  study  by  Rapp  (1982)  suggests  for  the  "best"  Airy-isostasy  a  compensation 
depth  of  50  km  rather  than  the  conventional  30  km,  a  probable  effect  of  the  trenches. 

The  primary  purpose  of  this  section  of  the  report  has  been  to  demonstrate 
the  practical  use  of  the  SYNBAPS-bathymetry  data  for  ocean  geoid  studies.  The 
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primitive  approach  taken  here  should  be  improved  for  future  studies,  and  of  par¬ 
ticular  interest  would  be  to  utilize  the  global  coverage  of  the  data  to  make  a 
comprehensive  analysis  of  global  oceanic  isostasy,  e.g.  through  the  use  of  empirical 
bathymetry  (geoid  transfer  functions),  class  fying  different  tectonic  ocean  areas. 
Results  of  such  an  analysis  could  then  be  used  for  the  important  "inverse"  problem  - 
to  determine  bathymetry  from  satellite  altimt?try  in  unsurveyed  areas. 

11.  Summary  and  Conclusions 

In  the  present  report,  many  different  topics  have  been  treated.  In  spite 
of  the  somewhat  diverse  composition,  I  hope  that  the  reader  still  has  felt  some 
kind  of  continuity  in  the  contents  -  the  idee  was  to  provide  the  general  background 
of  gravity  field  modelling  using  topographic/geologic  information,  stress  the 
similarities  to  and  possible  uses  of  geophysical  inversion  methods,  stress  the 
practical  benefits  of  using  spherical  harmonic  reference  fields  and  then  finally 
go  into  details  on  terrain  reductions  and  provide  the  theoretical  background  for 
the  FORTRAN  programme  published  with  this  report. 

In  the  first  part  of  the  report  (sections  2-5),  the  necessary  theoretical 
background  was  outlined,  including  the  concept  of  density  anomalies.  Opposed 
to  physical  densities,  density  anomalies  attain  both  positive  and  negative  values. 

In  spite  of  the  basic  ambiguity  of  potential  field  theory,  it  is  still  very  meaning¬ 
ful  to  work  with  density  anomalies,  even  without  having  defined  the  "normal"  density 
distribution  explicitly:  it  may  simply  be  taken  as  an  average  "expected"  geo¬ 
physical  model. 

Known  density  anomalies  -  topography,  isostasy,  geology  etc.  -  may  be  taken 
into  account  by  a  "remove-restore"  technique.  Then  values  of  anomalous  density 
etc.  are  assumed  to  be  known.  For  unknown  densities,  geophysical  inversion 
methods  may  be  used  to  provide  better  models  of  geologic  structures  or  to  provide 
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"optimal"  density  anomaly  values,  e.g.  finding  a  best  density  for  terrain  reduc¬ 
tions.  Furthermore,  it  was  outlined  how  the  "ultimate"  combination  of  generalized 
inversion  and  gravity  field  modelling  by  collocation  in  principle  could  be  done, 
to  allow  incorporation  of  e.g.  "non-gravity"  geophysical  information  in  gravity 
field  modelling.  In  the  second  part  of  the  report  (sections  6-7),  emphasis  was 
on  terrain  reductions.  Formulas  for  the  gravity  field  around  the  rectangular 
prism  -  the  natural  "building  stone"  when  elevations  are  given  as  digital  terrain 
models  -  were  given  and  evaluated.  The  .basic  "types"  of  terrain  reductions  were 
then  reviewed: 

-  Bouguer  reduction,  removal  of  the  visual  topography 

-  Isostatic  reductions,  removal  of  the  visual  topography  and  the  isostatic 
compensation 

-  Residual  reductions,  removal  only  of  the  short-wavelength,  "noisy"  topo¬ 
graphy 

The  terrain  correction,  frequently  treated  in  the  literature,  should  not  be 
viewed  as  a  terrain  reduction,  but  rather  just  as  an  important  mathematical 
auxiliary  quantity.  The  programme  "TC"  -  given  in  the  appendix  -  may  be  used 
to  compute  any  of  these  types  of  reductions  for  gravity,  deflections  and  height 
anomalies.  The  development  and  implementation  of  "TC"  represents  the  bulk  of 
"practical"  work  associated  with  this  report. 

For  error  studies  and  FFT-methods,  an  approximative  terrain  reduction  for¬ 
mulation  -  the  "linear  approximation"  -  plays  a  key  role.  The  accuracy  of  the 
approximation  was  investigated  both  for  a  theoretical  model  and  for  actual  data, 
with  the  conclusion  that  the  approximation  is  usually  acceptable.  However,  care 
should  be  exercised  in  alpine  areas,  especially  for  deflections  of  the  vertical. 

Assuming  the  validity  of  the  linear  approximation,  effective  FFT-methods 
for  terrain  effect  computations  were  outlined,  and  finally  a  key  topic  -  DTM 
resolution  requirements  -  were  studied.  Curves  were  given  to  compute  r.m.s. 
errors  for  terrain  effects  for  gravity,  deflections  and  height  anomalies,  based 
on  spacing  of  the  given  elevations  and  on  topography  covariance  functions. 


Finally  in  the  report,  empirical  data  were  investigated,  to  give  actual 
empirical  information  on 

-  topography  covariance  functions  for  various  types  of  areas  (they  turned  out 
to  be  generally  exponential) 

-  magnitudes  of  terrain  corrections 

-  magnitudes  of  RTM  geoid  effects  (they  are  actually  so  small  that  the 
indirect  effect  on  gravity  anomalies  may  be  neglected) 

-  degree  of  smoothing  obtained  using  terrain  reductions  for  actual  gravity 
data  (as  expected  variance  decreased,  correlation  length  increased  and 
anisotropy  was  diminished) 

-  relationship  of  altimetry  derived  ocean  geoids  to  geoids  computed  solely 
from  bathymetric  data  (including  bathymetry  covariance  functions). 

The  main  empirical  results  -  covariance  functions  for  topography  and  gravity  - 

are  contained  in  Figures  23-24,  33-34  and  Tables  4  and  6. 

Unfortunately  the  duration  of  the  author’s  stay  at  The  Ohio  State  University 
was  too  short  to  allow  the  implementation  and  practical  evaluation  of  some  key 
topics  of  the  report.  Therefore,  natural  extensions  of  the  present  study  would 
be 

-  implementation  of  hybrid  geophysical  inversion/collocation,  with  test  in 
an  area  of  well-known  geology  with  large  density  anomalies,  e.g.  a  salt- 
dome  province,  shelf  area,  etc. 

-  implementation  of  a  "hierarchial "  FFT  terrain  effect  computation  system, 
with  comparison  to  results  e.g.  obtained  by  "TC"  in  a  suitable  test  area. 

Other  directions  for  future  research  could  be  the  extension  of  the  analysis  of 

the  covariance  functions  for  topography  and  gravity,  to  include  more  regions  of 

different  types.  Ideally,  a  classification  of  "covariance  provinces",  based  e.g. 

on  existing  geographical  landscape  classification  systems,  could  be  attempted. 

The  availability  of  the  global  detailed  bathymetric  data  set  SYNBAPS  opens 
as  just  mentioned  possibilities  for  extensive  studies  of  the  relationship  between 
the  ocean  geoid  and  the  bathymetry.  Ocean  isostasy  could  be  studied  for  many 
types  of  areas  through  (FFT)  cross-covariance  analysis.  With  the  results  it 
could  then  e.g.  be  possible  to  design  optimum  filters  to  derive  the  bathymetry 
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from  satellite  altimetry  in  poorly  surveyed  areas.  Another  application  would 
be  to  "enhance"  geoid  variations  not  associated  directly  with  bathymetry,  such 
isostatic  geoid  anomalies  being  obviously  of  great  interest  for  geophysicists. 


APPENDIX 

C  -  A  Program  for  General  Terrain  Effect  Conmutations 


The  terrain  effect  computations  mentioned  in  the  present  report  have  been 
done  using  the  Fortran  Program  "TC",  outlined  in  the  sequel.  This  program  -  which 
is  written  in  structured  Fortran  (FORTRAN  77)  -  has  been  developed  and  implemented 
at  the  Amdahl  system  at  The  Ohio  State  University,  and  has  been  tested  against 
“synthetic"  topography  (cones,  cylinder  segments)  and  results  from  an  older,  some¬ 
what  different  ALGOL  program  at  the  Danish  Geodetic  Institute. 

The  program  uses  a  set  of  digital  terrain  models  (DTM's)  to  calculate  gravity, 
deflections  of  the  vertical  or  height  anomalies  at  specific  points,  using  the 
formulas  of  the  rectangular  prism,  outlined  in  section  6.1.  Four  different  mass 
model  types  may  be  specified:  topography,  topography  and  isostatic  compensation, 
terrain  corrections  or  residual  terrain  effects.  The  computations  may  either 
be  done  out  to  a  certain  distance  R,  or  a  specified  fixed  area  may  be  taken  into 
account.  A  fixed  density  of  2.67  for  the  topography  (h>0)  and  1.64  for  the  bathy¬ 
metry  (h<0)  is  used  presently,  but  it  may  be  changed  easily  in  the  start  of  the 
program.  For  the  isostatic. compensation,  an  Airy  isostasy  with  Moho  density  con¬ 
trast  0.4  g/cm3  and  normal  crust  thickness  32  km  is  used,  in  accordance  with  accepted 
"best"  continental  values. 

The  curvature  of  the  earth  is  taken  into  account  to  the  second  order,  through 
the  use  of  the  "super  elevation". 


SUPELV  =  H 


earth 


which  gives  the  z-shift  of  a  prism  in  distance  s  .below  the  tangential  plane  at 
the  computation  station  S.  This  approximation  is  valid  for  distances  up  to  several 
thousand  kilometers.  The  superelevation  is  also  utilized  for  terrain  corrections, 
so  for  very  large  computation  radii  terrain  corrections  must  be  combined  with 


Figure  48  The  use  of  the  two  DTM's  in  "TC"  terrain  effect  computations.  Elevations 
from  the  detailed  inner  grid  is  used  in  the  outlined  rectangle  covering  the  circle 
of  radius  Rx ,  centered  at  the  computation  point  S,  the  coarse  grid  being  utilized 
outside  this  rectangle.  In  the  inner-most  3x3  grid  elements  (the  "inner-zone") 
the  digital  terrain  model  may  be  densified  using  a  bicubic  spline  interpolation. 

the  spherical  Bouguer  plate  formulas  rather  than  just  the  plane  formula  ( 2it Gp h ) 
to  give  the  complete  topographic  effects. 

Two  digital  terrain  models  are  ordinarily  used  for  the  terrain  effect  computa¬ 
tions  in  "TC"  (although  both  are  not  necessary).  A  detailed  grid  (say  1  km  x  1  km 
point  elevations  )  is  used  out  to  a  computation  distance  Rj ,  and  then  a  coarse 
grid  (say  10x10  km  mean  elevations)  is  used  for  the  remaining  "remote  zones", 
see  Figure  48.  A  distance  specification  Rt  of  the  order  of  magnitude  twice 
the  outer  grid  spacing  is  usually  sufficient.  In  addition  to  the  "detailed"  and 
"coarse"  DTM,  an  additional  reference  DTM  must  be  specified  for  residual  terrain 
reductions.  This  DTM  is  used  together  with  a  parabolic  hyperboloid  (bilinear) 
interpolation  scheme  to  define  the  mean  elevation  surface. 


Each  DTM  must  be  stored  in  a  separate  file,  for  which  a  simple  standard  format 
is  used.  A  DTM  file  is  initiated  by  a 

Label:  <p2,  xx,  \2,  A4>,  aa  (REALM) 

defining  the  coverage  and  grid  spacing  of  the  DTM,  followed  by  the 

Elevations:  hlt  h2,  h3,  h4,  ...  ( INTEGER*2) 

• 

scanned  in  west  to  east  stripes,  starting  from  the  north.  Each  "stripe"  is  one 
record.  The  first  elevation  in  the  sequence  is  thus  the  NW-point,  the  last  the 
SE-point.  The  limits  of  the  DTM  (44,  <p2,  \lt  x2)  is  specified  by  geographical 
coordinates  of  the  outer  limits  of  the  rectangular  grid  "compartments"  (viewed 
as  mean  elevations),  see  Figure  47.  Thus  for  point  elevations,  half  the  grid 
spacing  must  be  added/subtracted  from  the  limits  to  get  the  correct  limits  for  the 
label . 

For  the  inner  grid,  an  alternative  NGS  format  may  be  used.  This  format, 
used  for  data  obtained  from  the  National  Geodetic  Survey,  consists  of  0.5'x0.5' 
point  elevations,  partitioned  in  7.5'xl°  blocks  each  containing  1800  elevations, 
the  blocks  stored  sequentially  as 

<p,  A  (REALM),  hlt  h2,  ...,  hl800  ( INTEGER*2) 
where  U,  x)  specifies  the  SE-point  of  the  block. 

The  program  demands  the  inner  and  outer  grids  to  be  consistent,  that  is, 
the  outer  grid  spacings  must  be  multipla  of  the  inner  grid  spacings,  and  the  gridlines 
of  the  outer  grid  must  be  grid  lines  also  in  the  inner  grid  (ideally  the  outer 
grid  should  be  simple  averages  of  the  inner  grid).  If  this  is  not  the  case  the 
program  will  terminate  with  an  error  message.  This  will  also  happen  if  too  small 
arrays  are  available  for  storing  elevations.  (The  program  only  stores  elevations 
for  the  smallest  possible  area  internally,  but  the  dynamic  storage  is  done  common 
for  all  stations  to  save  time  and  I/O  transfers.  It  will  therefore  be  advantageous 
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to  "block"  the  computations,  reducing  widely  separated  station  groups  in  different 
runs) . 

Each  grid  "compartment"  and  its  elevation  defines  a  rectangular  prism,  the 
effect  of  which  is  summed  up  to  give  the  complete  terrain  effect.  The  shift  between 
exact  and  approximative  prism  formulas  is  automatic,  defined  through  the  ratio 
distance  to  prism/grid  spacing,  cf.  Figure  7.  The  speed  of  the  program  may  be 
increased  by  reducing  this  ratio  ("R2EXAC",  "R2MACM"‘in  subroutine  "PRISMl"),  however 
at  a  price  of  degraded  computation  accuracy.  The  values  used  in  the  program  repre¬ 
sent  a  reasonable  trade-off,  determined  from  Figure  7  and  practical  experiments. 


S 


Figure  49  Spline  interpolation  of  elevations 
in  an  inner  zone  and  possible  modification 
to  give  the  "correct"  elevation  at  computation 
point  s. 


The  "innermost"  topography,  surrounding  the  computation  point,  is  of  critical 
importance  for  both  gravity  and  deflections.  Through  an  input  variable  "IZCODE" 
it  is  possible  to  govern  how  the  inner  zone  (3x3  elements,  cf.  Figure  48)  is  taken 
into  account.  For  stations  at  altitude  no  special  actions  will  usually  be  needed 
(IZCODE  3),  otherwise  the  elevations  are  interpolated  using  a  bicubic  spline  inter¬ 
polation  to  dense,  non-equidistant  inner-zone  grid,  used  for  summing  up  inner-zone 
effects.  When  a  station  is  known  to  be  on  the  earth's  surface  and  has  a  given 
elevation,  a  discrepancy  between  the  "model"  elevation  and  actual  elevation  at 
S  is  unavoidable.  For  deflections  and  undulations  this  discrepancy  is  unimportant, 
and  the  terrain  effect  should  be  computed  at  the  "model"  elevation  (IZCODE  0  or 
2).  For  gravity,  however,  the  topographic  or  RTM-effects  are  correlated  directly 


with  the  model  elevation,  so  in  this  case  the  terrain  model  should  be  modified 
to  reproduce  the  elevation  of  S  (IZCODE  1).  This  modification  is  done  smoothly 
within  the  inner-zone  (Figure  49),  assuming  the  discrepancy  to  be  due  to  erroneous 
(biased)  DTM  elevations.  Since  the  physical  mass  model  is  changed  independently 
for  each  individual  station,  some  care  should  be  exercised  when  using  this  option, 
especially  for  large  values  of  the  DTM  grid  spacings. 

In  the  sequel  the  programme  is  listed,  more  detailed  information  may  be  found 
in  the  introductionary  or  "en-route"  comments.  The  programme  is  modular  constructed, 
with  subroutines  "TC"  giving  effect  in  one  station,  "DTC"  effect  for  one  grid 
compartment  and  "PRISM  1"  giving  prism  formulas.  Station  input/output  should  be 
modified  (in  the  boginning  of  the  program)  to  satisfy  the  particular  needs  of 
the  user.  Typical  computation  times  at  the  OSU  Amdahl  system  is  around  4-5  points/ 
second.  An  input  example  is  given  at  the  end  of  the  program  listing. 
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